Алговики - Вклад [ru]

Участник:Роман Землянский/Zavolskov/Zemlyanskiy

2016-11-01T20:53:35Z

Bobru:

'''Метод «разделяй и властвуй» вычисления собственных значений и векторов симметричной трехдиагональной матрицы'''

Авторы статьи:

'''Завольсков Владислав'''

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

[[#Макроструктура алгоритма]]

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

[[#Информационный граф]]

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

'''Землянский Роман'''

[[#Общее описание алгоритма]]

[[#Математическое описание алгоритма]]

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

[[#Ресурс параллелизма алгоритма]]

[[#Входные и выходные данные алгоритма]]

[[#Свойства алгоритма]]

614 группа

= ЧАСТЬ. Свойства и структура алгоритмов =
== Общее описание алгоритма ==
'''Метод разделяй и властвуй вычисления собственных значений и векторов трёхдиагональной матрицы''' - это наиболее быстрый из существующих методов, если нужны все собственные значения и собственные векторы трехдиагональной матрицы, начиная с порядка n, примерно равного 26. (Точное значение этого порогового порядка зависит от компьютера.) Его численно устойчивая реализация весьма не тривиальна. В самом деле, хотя впервые метод был предложен еще в 1981 г., "правильный" способ его реализации был найден лишь в 1992 г. Этот способ воплощен LAPACK-программами ssyevd (для плотных матриц) и sstevd (для трехдиагональных матриц). В них стратегия "разделяй-и-влавствуй" используется для матриц порядка, большего чем 25. Для матриц меньшего порядка (или если нужны только собственные значения) происходит автоматический переход к QR-итерации.

== Математическое описание алгоритма ==
Пусть

:<math>
L = \begin{bmatrix}
a_{1} & b_{1}&&&&& \\
b_{1} & \ddots & \ddots \\
& \ddots & a_{m-1} & b_{m-1} \\
&& b_{m-1} & a_{m} & b_{m} \\
&&& b_{m} & a_{m+1} & b_{m+1} \\
&&&& b_{m+1} & \ddots \\
&&&&&& \ddots & b_{n-1} \\
&&&&&& b_{n-1} & a_{n} \\
\end{bmatrix}
= \begin{bmatrix}
a_{1} & b_{1} &&&&& \\
b_{1} & \ddots & \ddots \\
& \ddots & a_{m-1} & b_{m-1} \\
&& b_{m-1} & a_{m} - b_{m} \\
&&&& a_{m+1} - b_{m} & b_{m+1} \\
&&&& b_{m+1} & \ddots \\
&&&&&& \ddots & b_{n-1} \\
&&&&&& b_{n-1} & a_{n} \\
\end{bmatrix} +
</math>

<math>
+ \begin{bmatrix}
&&&&& \\
&&b_{m} & b_{m} \\
&&b_{m} & b_{m} \\
&&&&& \\
\end{bmatrix}

= \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m} * \begin{bmatrix} 0 \\ \vdots \\ 0 \\ 1 \\ 1 \\ 0 \\ \vdots \\ 0 \end{bmatrix} \begin{bmatrix} 0 , \ldots , 0 , 1 , 1 , 0 \ldots , 0 \end{bmatrix}

\equiv \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T}
</math>

Предположим, что нам известны спектральные разложения матриц <math>T_{1}</math> и <math> T_{2} </math>:
<math> T_{i} = Q_{i} \Lambda_{i} Q_{i}^{T} </math>. В действительности, они будут рекурсивно вычисляться тем же самым алгоритмом. Установим связь между собственными значениями матрицы Т и собственными значениями матриц
<math>T_{1}</math> и <math>T_{2}</math>. Имеем:

<math>
T = \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T}
= \begin{bmatrix} Q_{1} \Lambda_{1} Q_{1}^{T} & 0 \\ 0 & Q_{2} L_{2} Q_{2}^{T}\end{bmatrix} + b_{m}vv^{T}
= \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}(\begin{bmatrix} \Lambda_{1} & \\ & \Lambda_{2}\end{bmatrix} + b_{m}vv^{T})\begin{bmatrix} Q_{1}^{T} & 0 \\ 0 & Q_{2}^{T}\end{bmatrix}
</math>,

где

<math>
u = \begin{bmatrix} Q_{1}^{T} & 0 \\ 0 & Q_{2}^{T}\end{bmatrix}v
</math>

так как <math>v = \begin{bmatrix} 0 , \ldots , 0 , 1 , 1 , 0 \ldots , 0 \end{bmatrix}^T</math>, получим матрицу, состоящую из последнего столбца матрицы <math> Q_{1}^{T}</math> и первого столбца матрицы <math> Q_{2}^{T}</math>.

Следовательно, <math>T</math> имеет те же собственные значения, что и пдобная ей матрица <math>D + \rho uu^{T}</math>,
где <math>D = \begin{bmatrix} L_{1} & 0 \\ 0 & L_{2}\end{bmatrix}</math> - диагональная матрица,
<math>\rho = b_{m}</math> - число, а <math>u</math> - вектор.

Будем предполагать, не ограничивая общности, что диагональные элементы <math>d_{1}, \ldots, d_{n}</math> матрицы <math>D</math> упорядочены по убыванию: <math>d_{n} <= \ldots <=d_{1}</math>.

Чтобы найти собственные значения матрицы <math>D + \rho uu^{T}</math>, вычислим её характеристический многочлен, считая пока матрицу <math>D - \lambda I</math> невырожденной.
Тогда

<math>det(D + \rho uu^{T} - \lambda I) = det((D - \lambda I)(I + \rho (D- \lambda I)^{-1} uu^{T}))</math>.

Поскольку <math>D - \lambda I</math> невырожденна, <math>det(I + \rho (D - \lambda I)^{-1}uu^{T}) = 0</math> тогда и только тогда, когда <math>\lambda</math> - собственное значение. Заметим, что матрица <math>I + \rho (D - \lambda I)^{-1}uu^{T}</math> получается из единичной добавлением матрицы ранга 1. Определитель такой матрицы легко вычислить.

'''Лемма 1.''' Справедливо равенство <math>det(I + xy^{T}) = 1 + y^{T}x</math>, где <math>x</math> и <math>y</math> - векторы.

Таким образом,

<math>det(I + \rho (D - \lambda I)^{-1}uu^{T}) = 1 + \rho u^{T}(D - \lambda I)^{-1}u</math> <math> = 1 + \rho \sum_{i=1, n} \frac{u_{i}^{2}} {d_{i}-\lambda} \equiv f(\lambda)</math> ,

т.е. собственные значения матрицы <math>T</math> есть корни так называемого векового уравнения <math>f(\lambda) = 0</math>. Если все числа <math>d_{i}</math> различны и все <math>u_{i} <> 0</math> (случай общего положения), то <math>f(\lambda)</math> имеет график типа показанного на рис.1(где <math>n = 4</math> и <math>\rho > 0</math>).

[[File:Graphics.PNG|thumb|center|800px|Рис. 1. График функции <math> f(\lambda) = 1 + \frac{0.5}{1 - \lambda} + \frac{0.5}{2 - \lambda} + \frac{0.5}{3 - \lambda} + \frac{0.5}{4 - \lambda}</math>]]

Можно видеть, что прямая <math>y = 1</math> является горизонтальной асимптотой для этого графика, а прямые <math>\lambda = d_{i}</math> есть вертикальные асимптоты. Поскольку <math>f^{'}(\lambda) = \rho \sum_{i=1, n} \frac{u_{i}^{2}} {(d_{i}-\lambda)^{2}}> 0 </math>, функция возрастает всюду, кроме точек <math>\lambda = d_{i}</math>. Поэтому корни функции разделяются числами <math>d_{i}</math> и ещё один корень находится справа от точки <math>d_{1}</math> (на рис. 1 <math>d_{1} = 4</math>). (При <math>\rho<0</math> функция <math>f(\lambda)</math> всюду убывает и соответствующий корень находится слева от точки <math>d_{n}</math>). Для функции <math>f(\lambda)</math>, монотонной и гладкой на каждом из интервалов <math>(d_{i+1},d_{i})</math>, можно найти вариант метода Ньютона, который быстро и монотонно сходится к каждому из корней, если начальная точка взята в <math>(d_{i+1},d_{i})</math>. Нам достаточно знать, что на практике метод сходится к каждому собственному значению за ограниченное число шагов. Поскольку вычисление <math>f(\lambda)</math> и <math>f^{'}(\lambda)</math> стоит <math>O(n)</math> флопов, для вычисления одного собственного значения достаточно <math>O(n)</math>, а для вычисления всех <math>n</math> собственных значений матрицы <math>D + \rho uu^{T}</math> требуется <math>O(n^{2})</math> флопов.
Для собственных векторов матрицы <math>D + \rho uu^{T}</math> мы легко можем получить явные выражения.

'''Лемма 2.''' Если <math>\alpha</math> - собственное значение матрицы <math>D + \rho uu^{T}</math>, то соответствующий вектор равен <math>(D - \alpha I)^{-1}u</math>. Поскольку матрица <math>D - \alpha I</math> диагональная, для вычисления такого вектора достаточно <math>O(n)</math> флопов.

Доказательство.

<math>(D + \rho uu^{T})[(D - \alpha I)^{-1}u] = (D - \alpha I + \alpha I + \rho uu^{T})(D - \alpha I)^{-1}u</math>

<math>=u + \alpha (D - \alpha I)^{-1}u + u[\rho u^{T}(D - \alpha I)^{-1}u] </math>

<math>=u + \alpha(D - \alpha I)^{-1}u - u</math>

поскольку <math> \rho u^{T}(D - \alpha I)^{-1}u + 1 = f(\alpha) = 0 </math>

<math>=\alpha [(D - \alpha I)^{-1}u]</math>, что и требовалось.

Для вычисления по этой простой формуле всех <math>n</math> собственных векторов требуется <math>O(n^{2})</math> флопов. К сожалению, формула не обеспечивает численной устойчивости, так как для двух очень близких значений <math>\alpha_{i}</math> может давать неортогональные приближенные собственные векторы <math>u_{i}</math>. Потребовалось целое десятилетие для того, чтобы найти устойчивую альтернативу исходному описанию алгоритма. Снова детали будут обсуждаться позднее в данном разделе.

Алгоритм является рекурсивным.

=== Дефляция ===

До сих пор полагалось, что все <math>d_{i}</math> различны и все <math>u_{i}</math> отличны от нуля. Если это не так, вековое уравнение <math>f(\lambda)=0</math> имеет <math>k</math> вертикальных асимптот, где <math>k<n</math>, а потому <math>k</math> корней. Однако оказывается, что остальные <math>n - k</math> собственных значений могут быть определены без каких-либо усилий: если <math>d_{i}=d_{i+1}</math> или <math>u_{i}=0</math>, то легко показать, что <math>d_{i}</math> является собственным значением и для матрицы <math>D + \rho uu^{T}</math>. В такой ситуации мы говорим о ''дефляции''. На практике выбирается некоторое пороговое значение и дефляция для числа <math>d_{i}</math> регистрируется, если в смысле этого порога <math>d_{i}</math> достаточно близко к <math>d_{i+1}</math> либо <math>u_{i}</math> достаточно мало.

Основной выигрыш от использования дефляции состоит не в том, что убыстряется решение векового уравнения - этот этап в любом случае стоит лишь <math>O(n^{2})</math> операций. Выигрыш заключается в ускорении матричного умножения на последнем шаге алгоритма. Действительно, если <math>u_{i}=0</math>, то соответствующий собственный вектор есть i-й столбец <math>e_{i}</math> единичной матрицы. Это означает, что <math>e_{i}</math> является i-м столбцом в матрице <math>Q_{'}</math>, поэтому при формировании матрицы <math>Q</math> посредством левого умножения <math>Q_{1}</math> на <math>Q_{2}</math> вычисление i-го столбца не требует никаких затрат. Аналогичное упрощение имеет место в случае <math>d_{i} = d_{i+1}</math>. При дефляции многих собственных значений устраняется большая часть работы, связанной с матричным умножением.

== Вычислительное ядро алгоритма ==
Вычислительным ядром последовательной схемы решения является вычисление матрицы <math>Q</math> собственных векторов путём умножения матрицы <math>Q = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}</math> на матрицу <math>Q^{'}</math> Данная операция имеет сложность <math>cn^{3}</math> о чём и говорится в разделе [[#Последовательная сложность алгоритма]] .
Ей предшествует вычисление собственных значений и векторов матрицы <math> D + \rho uu^{T}</math>

== Макроструктура алгоритма ==
В разделе [[#Информационный граф]] описана структура алгоритма, в которой есть блок умножения матриц для вычисления собственных векторов, являющийся вычислительным ядром алгоритма. В соответствующем разделе ([[#Вычислительное ядро алгоритма]]) мы упоминали о том, что данному блоку предшествует вычисление собственных значений, которое производится [https://ru.wikipedia.org/wiki/%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%9D%D1%8C%D1%8E%D1%82%D0%BE%D0%BD%D0%B0 методом Ньютона].

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

Предположим, что некоторое <math>u_{i}</math>, хотя и мало, все же недостаточно мало для того, чтобы была зарегистрирована дефляция. В этом случае применение метода Ньютона к решению векового уравнения встречается с затруднениями. Вспомним, что пересчет приближённого решения <math>u_{j}</math> уравнения <math>f(\lambda) = 0</math> в методе Ньютона основан на следующих положениях:

1. Вблизи точки <math>\lambda = \lambda_{j}</math> функция <math>f(\lambda)</math> аппроксимируется линейной функцией <math>l(\lambda)</math>; график есть прямая линия, касающаяся графика функции <math>f(\lambda)</math> при <math>\lambda = \lambda_{j}</math>.

2. В качестве <math>\lambda_{j+1}</math> берётся нуль этого линейного приближения, т.е. <math>l(\lambda_{j+1})=0</math>.

Функция, показанная на рис.1, не доставляет видимых трудностей методу Ньютона, поскольку вблизи каждого своего нуля <math>f(\lambda)</math> достаточно хорошо аппроксимируется линейными функциями. Однако рассмотрим график функции на рис. 2. Она получена из функции на рис. 1 заменой значения .5 для <math>u_{i}^{2}</math> на .001. Это новое значение недостаточно мало для того, чтобы вызвать дефляцию. График функции в левой части рис.2 визуально не отличим от её вертикальных и горизонтальных асимптот, поэтому в правой части укрупненно воспроизведён фрагмент графика, прилегающий к вертикальной асимптоте <math>\lambda = 2</math>. Видно, что график слишком быстро "выполняет поворот" и для большей части значений <math>\lambda</math> почти горизонтален. Поэтому, применяя метод Ньютона почти к любому начальному приближению <math>\lambda_{0}</math>, мы получаем линейное приближение <math>l(\lambda)</math> с почти горизонтальным графиком и малым положительным угловым коэффициентом. В результате <math>\lambda_{1}</math> является отрицательным числом, огромным по абсолютной величине, которое совершенно бесполезно в качестве приближения к истинному корню.

[[File:GraphNewton.png|thumb|center|800px|Рис. 2. График функции <math> f(\lambda) = 1 + \frac{10^{-3}}{1 - \lambda} + \frac{10^{-3}}{2 - \lambda} + \frac{10^{-3}}{3 - \lambda} + \frac{10^{-3}}{4 - \lambda}</math>]]

Чтобы найти выход из этого положения, можно модифицировать метод Ньютона следующим образом: раз <math>f(\lambda)</math> нельзя хорошо приблизить линейной функцией <math>l(\lambda)</math>, попробуем взять в качестве приближения какую-нибудь другую простую функцию <math>h(\lambda)</math>. Нет ничего особого именно в прямых линиях: для метода Ньютона вместо <math>l(\lambda)</math> можно взять любое приближение <math>h(\lambda)</math>, значения и нули которого легко вычисляются. Функция <math>f(\lambda)</math> имеет полюсы в точках <math>d_{i}</math> и <math>d_{i+1}</math>, которые определяют её поведение в соответствующих окрестностях. Поэтому при поиске корня в интервале <math>(d_{i+1}, d_{i})</math> естественно выбрать функцию <math>h(\lambda)</math>, также имеющую эти полюсы, т.е. функцию вида
<math>h(\lambda)= \frac{c_{1}}{d_{i}-\lambda} + \frac{c_{2}}{d_{i+1}-\lambda} + c_{3}</math>

Константы <math>c_{1},c_{2}</math> и <math>c_{3}</math> обеспечивающие, что <math>h(\lambda)</math> есть приближение к <math>f(\lambda)</math>, можно выбрать несколькими способами. Отметим, что если <math>c_{1},c_{2}</math> и <math>c_{3}</math> уже известны, то уравнение <math>h(\lambda)=0</math> легко решается относительно <math>\lambda</math>, поскольку сводится к эквивалентному квадратному уравнению
<math>c_{1}(d_{i+1}-\lambda)+c_{2}(d_{i}-\lambda)+c_{3}(d_{i}-\lambda)(d_{i+1}-\lambda)=0</math>

Пусть <math>\lambda_{j}</math> - приближённое значение корня. определим <math>c_{1},c_{2}</math> и <math>c_{3}</math> так, чтобы
<math>\frac{c_{1}}{d_{i}-\lambda} + \frac{c_{2}}{d_{i+1}-\lambda} + c_{3} = h(\lambda) \approx f(\lambda) = 1 + \rho \sum_{k=1, n} \frac{u_{k}^{2}} {d_{k}-\lambda} </math>

для <math>\lambda</math> в окрестности <math>\lambda_{j}</math>. Заметим, что

<math>f(\lambda) = 1 + \rho \sum_{k=1, i} \frac{u_{k}^{2}} {d_{k}-\lambda} + \rho \sum_{k=i+1, n} \frac{u_{k}^{2}} {d_{k}-\lambda} \equiv 1 + \psi_{1}(\lambda) + \psi_{2}(\lambda)</math>.

Если <math>\lambda \in (d_{i+1},d_{i})</math>, то <math>\psi_{1}(\lambda)</math> есть сумма положительных слагаемых, а <math>\psi_{2}(\lambda)</math> - сумма отрицательных. Поэтому и <math>\psi_{1}(\lambda)</math>, и <math>\psi_{2}(\lambda)</math> могут быть вычислены с высокой точностью; однако при их сложении вполне вероятно взаимное уничтожение верных разрядов и потеря относительной точности в сумме. Возьмем числа <math>c_{1}</math> и <math>\hat{c_{1}}</math>, такие, что функция

<math>h_{1}(\lambda) \equiv \hat{c_{1}} + \frac{c_{1}}{d_{i}-\lambda}</math> удовлетворяет условиям <math>h_{1}(\lambda_{j}) = \psi_{1}(\lambda_{j})</math> и <math>h_{1}^{'}(\lambda_{j})=\psi_{1}^{'}(\lambda_{j})</math> (1)

Это означает, что гипербола, являющаяся графиком функции <math>h_{1}(\lambda)</math>, касается графика функции <math>\psi_{i}(\lambda)</math> при <math>\lambda = \lambda_{j}</math>. Два условия в (1) - это обычные условия метода Ньютона, за исключением того, что вместо прямой в качестве приближения используется гипербола. Легко проверить, что <math>c_{1}=\psi_{1}^{'}(\lambda_{j})(d_{i} - \lambda_{j})^{2}</math> и <math>\hat{c_{1}}=\psi_{1}(\lambda_{j}) - \psi_{1}^{'}(\lambda_{j})(d_{i} - \lambda_{j})</math>

Подобным же образом выбираем <math>c_{2}</math> и <math>\hat{c_{2}}</math> так, чтобы функция

<math>h_{2}(\lambda) \equiv \hat{c_{2}} + \frac{c_{2}}{d_{i+1}-\lambda}</math>

Удовлетворяла условиям

<math>h_{2}(\lambda_{j}) = \psi_{2}(\lambda_{j})</math> и <math>h_{2}^{'}(\lambda_{j})=\psi_{2}^{'}(\lambda_{j})</math> (2)

Наконец, полагаем

<math>h(\lambda) = 1 + h_{1}(\lambda) + h_{2}(\lambda) = (1 + \hat{c_{1}} + \hat{c_{2}} + \frac{c_{1}}{d_{i}-\lambda} + \frac{c_{2}}{d_{i+1}-\lambda} </math>

<math>\equiv c_{3} + \frac{c_{1}}{d_{i}-\lambda} + \frac{c_{2}}{d_{i+1}-\lambda}</math>.

== Схема реализации последовательного алгоритма ==
''Вычисление собственных значений и собственных векторов симметричной трехдиагональной матрицы посредством стратегии "разделяй и властвуй"'':

''proc dc_eig''<math>(T,Q,\Lambda)...</math> по входной матрице <math>T</math> вычисляются выходные матрицы <math>Q</math> и <math>\Lambda</math>, такие, что <math>T = Q\Lambda Q^{T}</math>

'''Если <math>T</math> - матрица размера <math>1</math> x <math>1</math>'''

1. Присвоить выходным параметрам значения <math> Q = 1, \Lambda = T</math>

'''Иначе'''

1. Представить <math>T</math> в виде <math> T = \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T} </math>

2. ''call dc_eig''<math>(T_{1},Q_{1},\Lambda_{1})</math>

3. ''call dc_eig''<math>(T_{2},Q_{2},\Lambda_{2})</math>

4. Построить <math>D+\rho uu^{T}</math> по <math> \Lambda_{1},\Lambda_{2}, Q_{1}, Q_{2}</math>

5. Найти матрицу собственных значений <math>\Lambda</math>

6. Найти матрицу собственных векторов <math>Q^{'}</math> для матрицы <math>D+\rho uu^{T}</math>

7. Построить матрицу собственных векторов <math>Q</math> для матрицы <math>T</math> : <math> Q = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}* Q^{'} </math>

8. Присвоить выходным параметрам значения <math>Q</math> и <math>\Lambda</math>

'''endif'''

== Последовательная сложность алгоритма ==
Проанализируем сложность алгоритма. Пусть <math>t(n)</math> - число флопов при обработке матрицы размера <math> n x n</math> процедурой ''dc_eig''. Тогда

<math> t(n) = 2t(n/2) </math> два рекурсивных обращения к ''dc_eig''<math>(T_{i},Q_{i},\Lambda_{i})</math>

<math> +O(n^{2})</math> вычисление собственных значений матрицы <math>D+\rho uu^{T}</math>

<math> +O(n^{2})</math> вычисление собственных векторов матрицы <math>D+\rho uu^{T}</math>

<math> +c*n^{3}</math> вычисление матрицы <math>Q = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}*Q^{'}</math>

Если <math> Q_{1}, Q_{2}</math> и <math>Q^{'}</math> рассматриваются как плотные матрицы и используется стандартный алгоритм матричного умножения, то константа <math> c </math> в последней строке равна 1. Таким образом, именно это умножение составляет наиболее трудоёмкую часть алгоритма в целом. Игнорируя члены порядка <math>n^{2}</math>, получаем <math>t(n) = 2t(n/2) + cn^{3}</math>. Решая это разностное уравнение, находим <math> t \approx c\frac{4}{3}n^{3} </math>

На практике константа <math>c</math> обычно гораздо меньше 1, потому что матрица <math>Q^{'}</math> весьма разрежена вследствие явления, называемого дефляцией.

== Информационный граф ==
В данном разделе представлен информационный граф алгоритма: на рисунке 3 изображена часть алгоритма, в которой происходит деление - "Divide", а на рисунке 4 описана часть "Conquer".

[[File:Divide.jpg|thumb|center|left|800px|Рис. 3. Информационный граф алгоритма "Разделяй и властвуй" - шаг Divide]]

[[File:ConquerPart.jpg|thumb|center|1200px|Рис. 4. Информационный граф алгоритма "Разделяй и властвуй" - шаг Conquer. Qk*k - матрица собственных значений размерности k. Lk*k - матрица собственных векторов размерности k.]]

Операции:

<math>D</math> - операция Divide разделения матрицы <math> T </math> на матрицы <math>T_{1}</math> и <math>T_{2}</math>

Индексы - размерность

<math>S_{k}</math> - операция Conquer сбора матрицы <math> D </math>

<math> L_{k} </math> - операция поиска собственных значений

<math> S^{'} </math> - операция поиска собственных векторов матрицы <math> Q^{'} </math>

<math> S </math> - операция поиска собственных векторов матрицы <math> Q </math>

== Ресурс параллелизма алгоритма ==
Рассмотрев информационный граф алгоритма, можно заметить, что в структуре имеется два параллельных блока - вызов рекурсивных функций ''dc_eig'' для вычисления собственных значений и векторов матриц <math>T_{1}</math> и <math>T_{2}</math> - это единственная часть алгоритма, в которой мы прибегаем к параллелизму. В случае реализации без параллелизма, функции ''dc_eig'' отрабатывают последовательно - сначала со входными параметрами <math>T_{1}, Q_{1}, \Lambda_{1}</math>, затем - <math>T_{2}, Q_{2}, \Lambda_{2}</math>.

== Входные и выходные данные алгоритма ==
На вход алгоритму подаётся трёхдиагональная матрица, описанная в разделе [[#Математическое описание алгоритма]]

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

== Свойства алгоритма ==
Описанный в данной статье алгоритм является самым быстрым алгоритмом, среди существующих: QR / Бисекция и обратная итерация / Разделяй и властвуй.

Последовательная сложность алгоритма: <math> t = c\frac{4}{3}n^{3} </math>

Параллельная сложность алгоритма (в силу использования дефляции): <math> O(n^{2.3)}</math> или в самых редких случаях <math> O(n^{2}) </math> (см. [3] стр. 8).

Особенности алгоритма:

1. Использование дефляции ([[#Дефляция]])

2. Использование адаптированного метода Ньютона ([[#Макроструктура алгоритма]])

= ЧАСТЬ. Программная реализация алгоритма =
Вторая часть описания алгоритмов в рамках AlgoWiki рассматривает все составные части процесса их реализации. Рассматривается как последовательная реализация алгоритма, так и параллельная. Описывается взаимосвязь свойств программ, реализующих алгоритм, и особенностей архитектуры компьютера, на которой они выполняются. Исследуется работа с памятью, локальность данных и вычислений, описывается масштабируемость и эффективность параллельных программ, производительность компьютеров, достигаемая на данной программе. Обсуждаются особенности реализации для разных классов архитектур компьютеров, приводятся ссылки на реализации в существующих библиотеках.

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

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

== Возможные способы и особенности параллельной реализации алгоритма ==
== Масштабируемость алгоритма и его реализации ==
TBD

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

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

== Существующие реализации алгоритма ==
Большого количества реализаций данного алгоритма найдено не было: единственный пример описан на страницах 7, 19 и 20 в [http://www.netlib.org/lapack/lawnspdf/lawn132.pdf этой статье], которая указана в списке литературы.

= Литература =
[1] Дж. Деммель, «Вычислительная линейная алгебра» //С. 230-235

[2] [http://www.netlib.org/utk/people/JackDongarra/PAPERS/104_1999_a-parallel-divide-and-conquer-algorithm.pdf Francoise Tisseury, Jack Dongarra, A Parallel Divide and Conquer Algorithm for the Symmetric Eigenvalue Problem on Distributed Memory Architectures]

[3] [http://www.netlib.org/lapack/lawnspdf/lawn132.pdf Francoise Tisseury, Jack Dongarra, Parallelizing the Divide and Conquer Algorithm for the Symmetric Tridiagonal Eigenvalue Problem on Distributed Memory Architectures]

[4] [https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_algorithm Алгоритм "Разделяй и властвуй" - Wikipedia]

[5] [http://www.cscamm.umd.edu/tadmor/pub/linear-stability/Gill_Tadmor_SISC90.pdf Doron Grill and Eitan Tadmor AN <math>O(N2)</math> METHOD FOR COMPUTING THE EIGENSYSTEM OF <math>N</math>x<math>N</math> SYMMETRIC TRIDIAGONAL MATRICES BY THE DIVIDE AND CONQUER APPROACH]

[[en:Description of algorithm properties and structure]]

Участник:Роман Землянский/Zavolskov/Zemlyanskiy

2016-11-01T20:48:33Z

Bobru:

'''Метод «разделяй и властвуй» вычисления собственных значений и векторов симметричной трехдиагональной матрицы'''

Авторы статьи:

'''Завольсков Владислав'''

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

[[#Макроструктура алгоритма]]

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

[[#Информационный граф]]

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

'''Землянский Роман'''

[[#Общее описание алгоритма]]

[[#Математическое описание алгоритма]]

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

[[#Ресурс параллелизма алгоритма]]

[[#Входные и выходные данные алгоритма]]

[[#Свойства алгоритма]]

614 группа

= ЧАСТЬ. Свойства и структура алгоритмов =
== Общее описание алгоритма ==
'''Метод разделяй и властвуй вычисления собственных значений и векторов трёхдиагональной матрицы''' - это наиболее быстрый из существующих методов, если нужны все собственные значения и собственные векторы трехдиагональной матрицы, начиная с порядка n, примерно равного 26. (Точное значение этого порогового порядка зависит от компьютера.) Его численно устойчивая реализация весьма не тривиальна. В самом деле, хотя впервые метод был предложен еще в 1981 г., "правильный" способ его реализации был найден лишь в 1992 г. Этот способ воплощен LAPACK-программами ssyevd (для плотных матриц) и sstevd (для трехдиагональных матриц). В них стратегия "разделяй-и-влавствуй" используется для матриц порядка, большего чем 25. Для матриц меньшего порядка (или если нужны только собственные значения) происходит автоматический переход к QR-итерации.

== Математическое описание алгоритма ==
Пусть

:<math>
L = \begin{bmatrix}
a_{1} & b_{1}&&&&& \\
b_{1} & \ddots & \ddots \\
& \ddots & a_{m-1} & b_{m-1} \\
&& b_{m-1} & a_{m} & b_{m} \\
&&& b_{m} & a_{m+1} & b_{m+1} \\
&&&& b_{m+1} & \ddots \\
&&&&&& \ddots & b_{n-1} \\
&&&&&& b_{n-1} & a_{n} \\
\end{bmatrix}
= \begin{bmatrix}
a_{1} & b_{1} &&&&& \\
b_{1} & \ddots & \ddots \\
& \ddots & a_{m-1} & b_{m-1} \\
&& b_{m-1} & a_{m} - b_{m} \\
&&&& a_{m+1} - b_{m} & b_{m+1} \\
&&&& b_{m+1} & \ddots \\
&&&&&& \ddots & b_{n-1} \\
&&&&&& b_{n-1} & a_{n} \\
\end{bmatrix} +
</math>

<math>
+ \begin{bmatrix}
&&&&& \\
&&b_{m} & b_{m} \\
&&b_{m} & b_{m} \\
&&&&& \\
\end{bmatrix}

= \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m} * \begin{bmatrix} 0 \\ \vdots \\ 0 \\ 1 \\ 1 \\ 0 \\ \vdots \\ 0 \end{bmatrix} \begin{bmatrix} 0 , \ldots , 0 , 1 , 1 , 0 \ldots , 0 \end{bmatrix}

\equiv \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T}
</math>

Предположим, что нам известны спектральные разложения матриц <math>T_{1}</math> и <math> T_{2} </math>:
<math> T_{i} = Q_{i} \Lambda_{i} Q_{i}^{T} </math>. В действительности, они будут рекурсивно вычисляться тем же самым алгоритмом. Установим связь между собственными значениями матрицы Т и собственными значениями матриц
<math>T_{1}</math> и <math>T_{2}</math>. Имеем:

<math>
T = \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T}
= \begin{bmatrix} Q_{1} \Lambda_{1} Q_{1}^{T} & 0 \\ 0 & Q_{2} L_{2} Q_{2}^{T}\end{bmatrix} + b_{m}vv^{T}
= \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}(\begin{bmatrix} \Lambda_{1} & \\ & \Lambda_{2}\end{bmatrix} + b_{m}vv^{T})\begin{bmatrix} Q_{1}^{T} & 0 \\ 0 & Q_{2}^{T}\end{bmatrix}
</math>,

где

<math>
u = \begin{bmatrix} Q_{1}^{T} & 0 \\ 0 & Q_{2}^{T}\end{bmatrix}v
</math>

так как <math>v = \begin{bmatrix} 0 , \ldots , 0 , 1 , 1 , 0 \ldots , 0 \end{bmatrix}^T</math>, получим матрицу, состоящую из последнего столбца матрицы <math> Q_{1}^{T}</math> и первого столбца матрицы <math> Q_{2}^{T}</math>.

Следовательно, <math>T</math> имеет те же собственные значения, что и пдобная ей матрица <math>D + \rho uu^{T}</math>,
где <math>D = \begin{bmatrix} L_{1} & 0 \\ 0 & L_{2}\end{bmatrix}</math> - диагональная матрица,
<math>\rho = b_{m}</math> - число, а <math>u</math> - вектор.

Будем предполагать, не ограничивая общности, что диагональные элементы <math>d_{1}, \ldots, d_{n}</math> матрицы <math>D</math> упорядочены по убыванию: <math>d_{n} <= \ldots <=d_{1}</math>.

Чтобы найти собственные значения матрицы <math>D + \rho uu^{T}</math>, вычислим её характеристический многочлен, считая пока матрицу <math>D - \lambda I</math> невырожденной.
Тогда

<math>det(D + \rho uu^{T} - \lambda I) = det((D - \lambda I)(I + \rho (D- \lambda I)^{-1} uu^{T}))</math>.

Поскольку <math>D - \lambda I</math> невырожденна, <math>det(I + \rho (D - \lambda I)^{-1}uu^{T}) = 0</math> тогда и только тогда, когда <math>\lambda</math> - собственное значение. Заметим, что матрица <math>I + \rho (D - \lambda I)^{-1}uu^{T}</math> получается из единичной добавлением матрицы ранга 1. Определитель такой матрицы легко вычислить.

'''Лемма 1.''' Справедливо равенство <math>det(I + xy^{T}) = 1 + y^{T}x</math>, где <math>x</math> и <math>y</math> - векторы.

Таким образом,

<math>det(I + \rho (D - \lambda I)^{-1}uu^{T}) = 1 + \rho u^{T}(D - \lambda I)^{-1}u</math> <math> = 1 + \rho \sum_{i=1, n} \frac{u_{i}^{2}} {d_{i}-\lambda} \equiv f(\lambda)</math> ,

т.е. собственные значения матрицы <math>T</math> есть корни так называемого векового уравнения <math>f(\lambda) = 0</math>. Если все числа <math>d_{i}</math> различны и все <math>u_{i} <> 0</math> (случай общего положения), то <math>f(\lambda)</math> имеет график типа показанного на рис.1(где <math>n = 4</math> и <math>\rho > 0</math>).

[[File:Graphics.PNG|thumb|center|800px|Рис. 1. График функции <math> f(\lambda) = 1 + \frac{0.5}{1 - \lambda} + \frac{0.5}{2 - \lambda} + \frac{0.5}{3 - \lambda} + \frac{0.5}{4 - \lambda}</math>]]

Можно видеть, что прямая <math>y = 1</math> является горизонтальной асимптотой для этого графика, а прямые <math>\lambda = d_{i}</math> есть вертикальные асимптоты. Поскольку <math>f^{'}(\lambda) = \rho \sum_{i=1, n} \frac{u_{i}^{2}} {(d_{i}-\lambda)^{2}}> 0 </math>, функция возрастает всюду, кроме точек <math>\lambda = d_{i}</math>. Поэтому корни функции разделяются числами <math>d_{i}</math> и ещё один корень находится справа от точки <math>d_{1}</math> (на рис. 1 <math>d_{1} = 4</math>). (При <math>\rho<0</math> функция <math>f(\lambda)</math> всюду убывает и соответствующий корень находится слева от точки <math>d_{n}</math>). Для функции <math>f(\lambda)</math>, монотонной и гладкой на каждом из интервалов <math>(d_{i+1},d_{i})</math>, можно найти вариант метода Ньютона, который быстро и монотонно сходится к каждому из корней, если начальная точка взята в <math>(d_{i+1},d_{i})</math>. Нам достаточно знать, что на практике метод сходится к каждому собственному значению за ограниченное число шагов. Поскольку вычисление <math>f(\lambda)</math> и <math>f^{'}(\lambda)</math> стоит <math>O(n)</math> флопов, для вычисления одного собственного значения достаточно <math>O(n)</math>, а для вычисления всех <math>n</math> собственных значений матрицы <math>D + \rho uu^{T}</math> требуется <math>O(n^{2})</math> флопов.
Для собственных векторов матрицы <math>D + \rho uu^{T}</math> мы легко можем получить явные выражения.

'''Лемма 2.''' Если <math>\alpha</math> - собственное значение матрицы <math>D + \rho uu^{T}</math>, то соответствующий вектор равен <math>(D - \alpha I)^{-1}u</math>. Поскольку матрица <math>D - \alpha I</math> диагональная, для вычисления такого вектора достаточно <math>O(n)</math> флопов.

Доказательство.

<math>(D + \rho uu^{T})[(D - \alpha I)^{-1}u] = (D - \alpha I + \alpha I + \rho uu^{T})(D - \alpha I)^{-1}u</math>

<math>=u + \alpha (D - \alpha I)^{-1}u + u[\rho u^{T}(D - \alpha I)^{-1}u] </math>

<math>=u + \alpha(D - \alpha I)^{-1}u - u</math>

поскольку <math> \rho u^{T}(D - \alpha I)^{-1}u + 1 = f(\alpha) = 0 </math>

<math>=\alpha [(D - \alpha I)^{-1}u]</math>, что и требовалось.

Для вычисления по этой простой формуле всех <math>n</math> собственных векторов требуется <math>O(n^{2})</math> флопов. К сожалению, формула не обеспечивает численной устойчивости, так как для двух очень близких значений <math>\alpha_{i}</math> может давать неортогональные приближенные собственные векторы <math>u_{i}</math>. Потребовалось целое десятилетие для того, чтобы найти устойчивую альтернативу исходному описанию алгоритма. Снова детали будут обсуждаться позднее в данном разделе.

Алгоритм является рекурсивным.

=== Дефляция ===

До сих пор полагалось, что все <math>d_{i}</math> различны и все <math>u_{i}</math> отличны от нуля. Если это не так, вековое уравнение <math>f(\lambda)=0</math> имеет <math>k</math> вертикальных асимптот, где <math>k<n</math>, а потому <math>k</math> корней. Однако оказывается, что остальные <math>n - k</math> собственных значений могут быть определены без каких-либо усилий: если <math>d_{i}=d_{i+1}</math> или <math>u_{i}=0</math>, то легко показать, что <math>d_{i}</math> является собственным значением и для матрицы <math>D + \rho uu^{T}</math>. В такой ситуации мы говорим о ''дефляции''. На практике выбирается некоторое пороговое значение и дефляция для числа <math>d_{i}</math> регистрируется, если в смысле этого порога <math>d_{i}</math> достаточно близко к <math>d_{i+1}</math> либо <math>u_{i}</math> достаточно мало.

Основной выигрыш от использования дефляции состоит не в том, что убыстряется решение векового уравнения - этот этап в любом случае стоит лишь <math>O(n^{2})</math> операций. Выигрыш заключается в ускорении матричного умножения на последнем шаге алгоритма. Действительно, если <math>u_{i}=0</math>, то соответствующий собственный вектор есть i-й столбец <math>e_{i}</math> единичной матрицы. Это означает, что <math>e_{i}</math> является i-м столбцом в матрице <math>Q_{'}</math>, поэтому при формировании матрицы <math>Q</math> посредством левого умножения <math>Q_{1}</math> на <math>Q_{2}</math> вычисление i-го столбца не требует никаких затрат. Аналогичное упрощение имеет место в случае <math>d_{i} = d_{i+1}</math>. При дефляции многих собственных значений устраняется большая часть работы, связанной с матричным умножением.

== Вычислительное ядро алгоритма ==
Вычислительным ядром последовательной схемы решения является вычисление матрицы <math>Q</math> собственных векторов путём умножения матрицы <math>Q = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}</math> на матрицу <math>Q^{'}</math> Данная операция имеет сложность <math>cn^{3}</math> о чём и говорится в разделе [[#Последовательная сложность алгоритма]] .
Ей предшествует вычисление собственных значений и векторов матрицы <math> D + \rho uu^{T}</math>

== Макроструктура алгоритма ==
В разделе [[#Информационный граф]] описана структура алгоритма, в которой есть блок умножения матриц для вычисления собственных векторов, являющийся вычислительным ядром алгоритма. В соответствующем разделе ([[#Вычислительное ядро алгоритма]]) мы упоминали о том, что данному блоку предшествует вычисление собственных значений, которое производится [https://ru.wikipedia.org/wiki/%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%9D%D1%8C%D1%8E%D1%82%D0%BE%D0%BD%D0%B0 методом Ньютона].

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

Предположим, что некоторое <math>u_{i}</math>, хотя и мало, все же недостаточно мало для того, чтобы была зарегистрирована дефляция. В этом случае применение метода Ньютона к решению векового уравнения встречается с затруднениями. Вспомним, что пересчет приближённого решения <math>u_{j}</math> уравнения <math>f(\lambda) = 0</math> в методе Ньютона основан на следующих положениях:

1. Вблизи точки <math>\lambda = \lambda_{j}</math> функция <math>f(\lambda)</math> аппроксимируется линейной функцией <math>l(\lambda)</math>; график есть прямая линия, касающаяся графика функции <math>f(\lambda)</math> при <math>\lambda = \lambda_{j}</math>.

2. В качестве <math>\lambda_{j+1}</math> берётся нуль этого линейного приближения, т.е. <math>l(\lambda_{j+1})=0</math>.

Функция, показанная на рис.1, не доставляет видимых трудностей методу Ньютона, поскольку вблизи каждого своего нуля <math>f(\lambda)</math> достаточно хорошо аппроксимируется линейными функциями. Однако рассмотрим график функции на рис. 2. Она получена из функции на рис. 1 заменой значения .5 для <math>u_{i}^{2}</math> на .001. Это новое значение недостаточно мало для того, чтобы вызвать дефляцию. График функции в левой части рис.2 визуально не отличим от её вертикальных и горизонтальных асимптот, поэтому в правой части укрупненно воспроизведён фрагмент графика, прилегающий к вертикальной асимптоте <math>\lambda = 2</math>. Видно, что график слишком быстро "выполняет поворот" и для большей части значений <math>\lambda</math> почти горизонтален. Поэтому, применяя метод Ньютона почти к любому начальному приближению <math>\lambda_{0}</math>, мы получаем линейное приближение <math>l(\lambda)</math> с почти горизонтальным графиком и малым положительным угловым коэффициентом. В результате <math>\lambda_{1}</math> является отрицательным числом, огромным по абсолютной величине, которое совершенно бесполезно в качестве приближения к истинному корню.

[[File:GraphNewton.png|thumb|center|800px|Рис. 2. График функции <math> f(\lambda) = 1 + \frac{10^{-3}}{1 - \lambda} + \frac{10^{-3}}{2 - \lambda} + \frac{10^{-3}}{3 - \lambda} + \frac{10^{-3}}{4 - \lambda}</math>]]

Чтобы найти выход из этого положения, можно модифицировать метод Ньютона следующим образом: раз <math>f(\lambda)</math> нельзя хорошо приблизить линейной функцией <math>l(\lambda)</math>, попробуем взять в качестве приближения какую-нибудь другую простую функцию <math>h(\lambda)</math>. Нет ничего особого именно в прямых линиях: для метода Ньютона вместо <math>l(\lambda)</math> можно взять любое приближение <math>h(\lambda)</math>, значения и нули которого легко вычисляются. Функция <math>f(\lambda)</math> имеет полюсы в точках <math>d_{i}</math> и <math>d_{i+1}</math>, которые определяют её поведение в соответствующих окрестностях. Поэтому при поиске корня в интервале <math>(d_{i+1}, d_{i})</math> естественно выбрать функцию <math>h(\lambda)</math>, также имеющую эти полюсы, т.е. функцию вида
<math>h(\lambda)= \frac{c_{1}}{d_{i}-\lambda} + \frac{c_{2}}{d_{i+1}-\lambda} + c_{3}</math>

Константы <math>c_{1},c_{2}</math> и <math>c_{3}</math> обеспечивающие, что <math>h(\lambda)</math> есть приближение к <math>f(\lambda)</math>, можно выбрать несколькими способами. Отметим, что если <math>c_{1},c_{2}</math> и <math>c_{3}</math> уже известны, то уравнение <math>h(\lambda)=0</math> легко решается относительно <math>\lambda</math>, поскольку сводится к эквивалентному квадратному уравнению
<math>c_{1}(d_{i+1}-\lambda)+c_{2}(d_{i}-\lambda)+c_{3}(d_{i}-\lambda)(d_{i+1}-\lambda)=0</math>

Пусть <math>\lambda_{j}</math> - приближённое значение корня. определим <math>c_{1},c_{2}</math> и <math>c_{3}</math> так, чтобы
<math>\frac{c_{1}}{d_{i}-\lambda} + \frac{c_{2}}{d_{i+1}-\lambda} + c_{3} = h(\lambda) \approx f(\lambda) = 1 + \rho \sum_{k=1, n} \frac{u_{k}^{2}} {d_{k}-\lambda} </math>

для <math>\lambda</math> в окрестности <math>\lambda_{j}</math>. Заметим, что

<math>f(\lambda) = 1 + \rho \sum_{k=1, i} \frac{u_{k}^{2}} {d_{k}-\lambda} + \rho \sum_{k=i+1, n} \frac{u_{k}^{2}} {d_{k}-\lambda} \equiv 1 + \psi_{1}(\lambda) + \psi_{2}(\lambda)</math>.

Если <math>\lambda \in (d_{i+1},d_{i})</math>, то <math>\psi_{1}(\lambda)</math> есть сумма положительных слагаемых, а <math>\psi_{2}(\lambda)</math> - сумма отрицательных. Поэтому и <math>\psi_{1}(\lambda)</math>, и <math>\psi_{2}(\lambda)</math> могут быть вычислены с высокой точностью; однако при их сложении вполне вероятно взаимное уничтожение верных разрядов и потеря относительной точности в сумме. Возьмем числа <math>c_{1}</math> и <math>\hat{c_{1}}</math>, такие, что функция

<math>h_{1}(\lambda) \equiv \hat{c_{1}} + \frac{c_{1}}{d_{i}-\lambda}</math> удовлетворяет условиям <math>h_{1}(\lambda_{j}) = \psi_{1}(\lambda_{j})</math> и <math>h_{1}^{'}(\lambda_{j})=\psi_{1}^{'}(\lambda_{j})</math> (1)

Это означает, что гипербола, являющаяся графиком функции <math>h_{1}(\lambda)</math>, касается графика функции <math>\psi_{i}(\lambda)</math> при <math>\lambda = \lambda_{j}</math>. Два условия в (1) - это обычные условия метода Ньютона, за исключением того, что вместо прямой в качестве приближения используется гипербола. Легко проверить, что <math>c_{1}=\psi_{1}^{'}(\lambda_{j})(d_{i} - \lambda_{j})^{2}</math> и <math>\hat{c_{1}}=\psi_{1}(\lambda_{j}) - \psi_{1}^{'}(\lambda_{j})(d_{i} - \lambda_{j})</math>

Подобным же образом выбираем <math>c_{2}</math> и <math>\hat{c_{2}}</math> так, чтобы функция

<math>h_{2}(\lambda) \equiv \hat{c_{2}} + \frac{c_{2}}{d_{i+1}-\lambda}</math>

Удовлетворяла условиям

<math>h_{2}(\lambda_{j}) = \psi_{2}(\lambda_{j})</math> и <math>h_{2}^{'}(\lambda_{j})=\psi_{2}^{'}(\lambda_{j})</math> (2)

Наконец, полагаем

<math>h(\lambda) = 1 + h_{1}(\lambda) + h_{2}(\lambda) = (1 + \hat{c_{1}} + \hat{c_{2}} + \frac{c_{1}}{d_{i}-\lambda} + \frac{c_{2}}{d_{i+1}-\lambda} </math>

<math>\equiv c_{3} + \frac{c_{1}}{d_{i}-\lambda} + \frac{c_{2}}{d_{i+1}-\lambda}</math>.

== Схема реализации последовательного алгоритма ==
''Вычисление собственных значений и собственных векторов симметричной трехдиагональной матрицы посредством стратегии "разделяй и властвуй"'':

''proc dc_eig''<math>(T,Q,\Lambda)...</math> по входной матрице <math>T</math> вычисляются выходные матрицы <math>Q</math> и <math>\Lambda</math>, такие, что <math>T = Q\Lambda Q^{T}</math>

'''Если <math>T</math> - матрица размера <math>1</math> x <math>1</math>'''

1. Присвоить выходным параметрам значения <math> Q = 1, \Lambda = T</math>

'''Иначе'''

1. Представить <math>T</math> в виде <math> T = \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T} </math>

2. ''call dc_eig''<math>(T_{1},Q_{1},\Lambda_{1})</math>

3. ''call dc_eig''<math>(T_{2},Q_{2},\Lambda_{2})</math>

4. Построить <math>D+\rho uu^{T}</math> по <math> \Lambda_{1},\Lambda_{2}, Q_{1}, Q_{2}</math>

5. Найти матрицу собственных значений <math>\Lambda</math>

6. Найти матрицу собственных векторов <math>Q^{'}</math> для матрицы <math>D+\rho uu^{T}</math>

7. Построить матрицу собственных векторов <math>Q</math> для матрицы <math>T</math> : <math> Q = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}* Q^{'} </math>

8. Присвоить выходным параметрам значения <math>Q</math> и <math>\Lambda</math>

'''endif'''

== Последовательная сложность алгоритма ==
Проанализируем сложность алгоритма. Пусть <math>t(n)</math> - число флопов при обработке матрицы размера <math> n x n</math> процедурой ''dc_eig''. Тогда

<math> t(n) = 2t(n/2) </math> два рекурсивных обращения к ''dc_eig''<math>(T_{i},Q_{i},\Lambda_{i})</math>

<math> +O(n^{2})</math> вычисление собственных значений матрицы <math>D+\rho uu^{T}</math>

<math> +O(n^{2})</math> вычисление собственных векторов матрицы <math>D+\rho uu^{T}</math>

<math> +c*n^{3}</math> вычисление матрицы <math>Q = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}*Q^{'}</math>

Если <math> Q_{1}, Q_{2}</math> и <math>Q^{'}</math> рассматриваются как плотные матрицы и используется стандартный алгоритм матричного умножения, то константа <math> c </math> в последней строке равна 1. Таким образом, именно это умножение составляет наиболее трудоёмкую часть алгоритма в целом. Игнорируя члены порядка <math>n^{2}</math>, получаем <math>t(n) = 2t(n/2) + cn^{3}</math>. Решая это разностное уравнение, находим <math> t \approx c\frac{4}{3}n^{3} </math>

На практике константа <math>c</math> обычно гораздо меньше 1, потому что матрица <math>Q^{'}</math> весьма разрежена вследствие явления, называемого дефляцией.

== Информационный граф ==
В данном разделе представлен информационный граф алгоритма: на рисунке 3 изображена часть алгоритма, в которой происходит деление - "Divide", а на рисунке 4 описана часть "Conquer".

[[File:Divide.jpg|thumb|center|left|800px|Рис. 3. Информационный граф алгоритма "Разделяй и властвуй" - шаг Divide]]

[[File:ConquerPart.jpg|thumb|center|1200px|Рис. 4. Информационный граф алгоритма "Разделяй и властвуй" - шаг Conquer. Qk*k - матрица собственных значений размерности k. Lk*k - матрица собственных векторов размерности k.]]

Операции:

<math>D</math> - операция Divide разделения матрицы <math> T </math> на матрицы <math>T_{1}</math> и <math>T_{2}</math>

Индексы - размерность

<math>S_{k}</math> - операция Conquer сбора матрицы <math> D </math>

<math> L </math> - операция поиска собственных значений

<math> S^{'} </math> - операция поиска собственных векторов матрицы <math> Q^{'} </math>

<math> S </math> - операция поиска собственных векторов матрицы <math> Q </math>

== Ресурс параллелизма алгоритма ==
Рассмотрев информационный граф алгоритма, можно заметить, что в структуре имеется два параллельных блока - вызов рекурсивных функций ''dc_eig'' для вычисления собственных значений и векторов матриц <math>T_{1}</math> и <math>T_{2}</math> - это единственная часть алгоритма, в которой мы прибегаем к параллелизму. В случае реализации без параллелизма, функции ''dc_eig'' отрабатывают последовательно - сначала со входными параметрами <math>T_{1}, Q_{1}, \Lambda_{1}</math>, затем - <math>T_{2}, Q_{2}, \Lambda_{2}</math>.

== Входные и выходные данные алгоритма ==
На вход алгоритму подаётся трёхдиагональная матрица, описанная в разделе [[#Математическое описание алгоритма]]

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

== Свойства алгоритма ==
Описанный в данной статье алгоритм является самым быстрым алгоритмом, среди существующих: QR / Бисекция и обратная итерация / Разделяй и властвуй.

Последовательная сложность алгоритма: <math> t = c\frac{4}{3}n^{3} </math>

Параллельная сложность алгоритма (в силу использования дефляции): <math> O(n^{2.3)}</math> или в самых редких случаях <math> O(n^{2}) </math> (см. [3] стр. 8).

Особенности алгоритма:

1. Использование дефляции ([[#Дефляция]])

2. Использование адаптированного метода Ньютона ([[#Макроструктура алгоритма]])

= ЧАСТЬ. Программная реализация алгоритма =
Вторая часть описания алгоритмов в рамках AlgoWiki рассматривает все составные части процесса их реализации. Рассматривается как последовательная реализация алгоритма, так и параллельная. Описывается взаимосвязь свойств программ, реализующих алгоритм, и особенностей архитектуры компьютера, на которой они выполняются. Исследуется работа с памятью, локальность данных и вычислений, описывается масштабируемость и эффективность параллельных программ, производительность компьютеров, достигаемая на данной программе. Обсуждаются особенности реализации для разных классов архитектур компьютеров, приводятся ссылки на реализации в существующих библиотеках.

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

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

== Возможные способы и особенности параллельной реализации алгоритма ==
== Масштабируемость алгоритма и его реализации ==
TBD

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

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

== Существующие реализации алгоритма ==
Большого количества реализаций данного алгоритма найдено не было: единственный пример описан на страницах 7, 19 и 20 в [http://www.netlib.org/lapack/lawnspdf/lawn132.pdf этой статье], которая указана в списке литературы.

= Литература =
[1] Дж. Деммель, «Вычислительная линейная алгебра» //С. 230-235

[2] [http://www.netlib.org/utk/people/JackDongarra/PAPERS/104_1999_a-parallel-divide-and-conquer-algorithm.pdf Francoise Tisseury, Jack Dongarra, A Parallel Divide and Conquer Algorithm for the Symmetric Eigenvalue Problem on Distributed Memory Architectures]

[3] [http://www.netlib.org/lapack/lawnspdf/lawn132.pdf Francoise Tisseury, Jack Dongarra, Parallelizing the Divide and Conquer Algorithm for the Symmetric Tridiagonal Eigenvalue Problem on Distributed Memory Architectures]

[4] [https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_algorithm Алгоритм "Разделяй и властвуй" - Wikipedia]

[5] [http://www.cscamm.umd.edu/tadmor/pub/linear-stability/Gill_Tadmor_SISC90.pdf Doron Grill and Eitan Tadmor AN <math>O(N2)</math> METHOD FOR COMPUTING THE EIGENSYSTEM OF <math>N</math>x<math>N</math> SYMMETRIC TRIDIAGONAL MATRICES BY THE DIVIDE AND CONQUER APPROACH]

[[en:Description of algorithm properties and structure]]

Участник:Роман Землянский/Zavolskov/Zemlyanskiy

2016-11-01T20:44:32Z

Bobru:

'''Метод «разделяй и властвуй» вычисления собственных значений и векторов симметричной трехдиагональной матрицы'''

Авторы статьи:

'''Завольсков Владислав'''

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

[[#Макроструктура алгоритма]]

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

[[#Информационный граф]]

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

'''Землянский Роман'''

[[#Общее описание алгоритма]]

[[#Математическое описание алгоритма]]

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

[[#Ресурс параллелизма алгоритма]]

[[#Входные и выходные данные алгоритма]]

[[#Свойства алгоритма]]

614 группа

= ЧАСТЬ. Свойства и структура алгоритмов =
== Общее описание алгоритма ==
'''Метод разделяй и властвуй вычисления собственных значений и векторов трёхдиагональной матрицы''' - это наиболее быстрый из существующих методов, если нужны все собственные значения и собственные векторы трехдиагональной матрицы, начиная с порядка n, примерно равного 26. (Точное значение этого порогового порядка зависит от компьютера.) Его численно устойчивая реализация весьма не тривиальна. В самом деле, хотя впервые метод был предложен еще в 1981 г., "правильный" способ его реализации был найден лишь в 1992 г. Этот способ воплощен LAPACK-программами ssyevd (для плотных матриц) и sstevd (для трехдиагональных матриц). В них стратегия "разделяй-и-влавствуй" используется для матриц порядка, большего чем 25. Для матриц меньшего порядка (или если нужны только собственные значения) происходит автоматический переход к QR-итерации.

== Математическое описание алгоритма ==
Пусть

:<math>
L = \begin{bmatrix}
a_{1} & b_{1}&&&&& \\
b_{1} & \ddots & \ddots \\
& \ddots & a_{m-1} & b_{m-1} \\
&& b_{m-1} & a_{m} & b_{m} \\
&&& b_{m} & a_{m+1} & b_{m+1} \\
&&&& b_{m+1} & \ddots \\
&&&&&& \ddots & b_{n-1} \\
&&&&&& b_{n-1} & a_{n} \\
\end{bmatrix}
= \begin{bmatrix}
a_{1} & b_{1} &&&&& \\
b_{1} & \ddots & \ddots \\
& \ddots & a_{m-1} & b_{m-1} \\
&& b_{m-1} & a_{m} - b_{m} \\
&&&& a_{m+1} - b_{m} & b_{m+1} \\
&&&& b_{m+1} & \ddots \\
&&&&&& \ddots & b_{n-1} \\
&&&&&& b_{n-1} & a_{n} \\
\end{bmatrix} +
</math>

<math>
+ \begin{bmatrix}
&&&&& \\
&&b_{m} & b_{m} \\
&&b_{m} & b_{m} \\
&&&&& \\
\end{bmatrix}

= \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m} * \begin{bmatrix} 0 \\ \vdots \\ 0 \\ 1 \\ 1 \\ 0 \\ \vdots \\ 0 \end{bmatrix} \begin{bmatrix} 0 , \ldots , 0 , 1 , 1 , 0 \ldots , 0 \end{bmatrix}

\equiv \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T}
</math>

Предположим, что нам известны спектральные разложения матриц <math>T_{1}</math> и <math> T_{2} </math>:
<math> T_{i} = Q_{i} \Lambda_{i} Q_{i}^{T} </math>. В действительности, они будут рекурсивно вычисляться тем же самым алгоритмом. Установим связь между собственными значениями матрицы Т и собственными значениями матриц
<math>T_{1}</math> и <math>T_{2}</math>. Имеем:

<math>
T = \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T}
= \begin{bmatrix} Q_{1} \Lambda_{1} Q_{1}^{T} & 0 \\ 0 & Q_{2} L_{2} Q_{2}^{T}\end{bmatrix} + b_{m}vv^{T}
= \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}(\begin{bmatrix} \Lambda_{1} & \\ & \Lambda_{2}\end{bmatrix} + b_{m}vv^{T})\begin{bmatrix} Q_{1}^{T} & 0 \\ 0 & Q_{2}^{T}\end{bmatrix}
</math>,

где

<math>
u = \begin{bmatrix} Q_{1}^{T} & 0 \\ 0 & Q_{2}^{T}\end{bmatrix}v
</math>

так как <math>v = \begin{bmatrix} 0 , \ldots , 0 , 1 , 1 , 0 \ldots , 0 \end{bmatrix}^T</math>, получим матрицу, состоящую из последнего столбца матрицы <math> Q_{1}^{T}</math> и первого столбца матрицы <math> Q_{2}^{T}</math>.

Следовательно, <math>T</math> имеет те же собственные значения, что и пдобная ей матрица <math>D + \rho uu^{T}</math>,
где <math>D = \begin{bmatrix} L_{1} & 0 \\ 0 & L_{2}\end{bmatrix}</math> - диагональная матрица,
<math>\rho = b_{m}</math> - число, а <math>u</math> - вектор.

Будем предполагать, не ограничивая общности, что диагональные элементы <math>d_{1}, \ldots, d_{n}</math> матрицы <math>D</math> упорядочены по убыванию: <math>d_{n} <= \ldots <=d_{1}</math>.

Чтобы найти собственные значения матрицы <math>D + \rho uu^{T}</math>, вычислим её характеристический многочлен, считая пока матрицу <math>D - \lambda I</math> невырожденной.
Тогда

<math>det(D + \rho uu^{T} - \lambda I) = det((D - \lambda I)(I + \rho (D- \lambda I)^{-1} uu^{T}))</math>.

Поскольку <math>D - \lambda I</math> невырожденна, <math>det(I + \rho (D - \lambda I)^{-1}uu^{T}) = 0</math> тогда и только тогда, когда <math>\lambda</math> - собственное значение. Заметим, что матрица <math>I + \rho (D - \lambda I)^{-1}uu^{T}</math> получается из единичной добавлением матрицы ранга 1. Определитель такой матрицы легко вычислить.

'''Лемма 1.''' Справедливо равенство <math>det(I + xy^{T}) = 1 + y^{T}x</math>, где <math>x</math> и <math>y</math> - векторы.

Таким образом,

<math>det(I + \rho (D - \lambda I)^{-1}uu^{T}) = 1 + \rho u^{T}(D - \lambda I)^{-1}u</math> <math> = 1 + \rho \sum_{i=1, n} \frac{u_{i}^{2}} {d_{i}-\lambda} \equiv f(\lambda)</math> ,

т.е. собственные значения матрицы <math>T</math> есть корни так называемого векового уравнения <math>f(\lambda) = 0</math>. Если все числа <math>d_{i}</math> различны и все <math>u_{i} <> 0</math> (случай общего положения), то <math>f(\lambda)</math> имеет график типа показанного на рис.1(где <math>n = 4</math> и <math>\rho > 0</math>).

[[File:Graphics.PNG|thumb|center|800px|Рис. 1. График функции <math> f(\lambda) = 1 + \frac{0.5}{1 - \lambda} + \frac{0.5}{2 - \lambda} + \frac{0.5}{3 - \lambda} + \frac{0.5}{4 - \lambda}</math>]]

Можно видеть, что прямая <math>y = 1</math> является горизонтальной асимптотой для этого графика, а прямые <math>\lambda = d_{i}</math> есть вертикальные асимптоты. Поскольку <math>f^{'}(\lambda) = \rho \sum_{i=1, n} \frac{u_{i}^{2}} {(d_{i}-\lambda)^{2}}> 0 </math>, функция возрастает всюду, кроме точек <math>\lambda = d_{i}</math>. Поэтому корни функции разделяются числами <math>d_{i}</math> и ещё один корень находится справа от точки <math>d_{1}</math> (на рис. 1 <math>d_{1} = 4</math>). (При <math>\rho<0</math> функция <math>f(\lambda)</math> всюду убывает и соответствующий корень находится слева от точки <math>d_{n}</math>). Для функции <math>f(\lambda)</math>, монотонной и гладкой на каждом из интервалов <math>(d_{i+1},d_{i})</math>, можно найти вариант метода Ньютона, который быстро и монотонно сходится к каждому из корней, если начальная точка взята в <math>(d_{i+1},d_{i})</math>. Нам достаточно знать, что на практике метод сходится к каждому собственному значению за ограниченное число шагов. Поскольку вычисление <math>f(\lambda)</math> и <math>f^{'}(\lambda)</math> стоит <math>O(n)</math> флопов, для вычисления одного собственного значения достаточно <math>O(n)</math>, а для вычисления всех <math>n</math> собственных значений матрицы <math>D + \rho uu^{T}</math> требуется <math>O(n^{2})</math> флопов.
Для собственных векторов матрицы <math>D + \rho uu^{T}</math> мы легко можем получить явные выражения.

'''Лемма 2.''' Если <math>\alpha</math> - собственное значение матрицы <math>D + \rho uu^{T}</math>, то соответствующий вектор равен <math>(D - \alpha I)^{-1}u</math>. Поскольку матрица <math>D - \alpha I</math> диагональная, для вычисления такого вектора достаточно <math>O(n)</math> флопов.

Доказательство.

<math>(D + \rho uu^{T})[(D - \alpha I)^{-1}u] = (D - \alpha I + \alpha I + \rho uu^{T})(D - \alpha I)^{-1}u</math>

<math>=u + \alpha (D - \alpha I)^{-1}u + u[\rho u^{T}(D - \alpha I)^{-1}u] </math>

<math>=u + \alpha(D - \alpha I)^{-1}u - u</math>

поскольку <math> \rho u^{T}(D - \alpha I)^{-1}u + 1 = f(\alpha) = 0 </math>

<math>=\alpha [(D - \alpha I)^{-1}u]</math>, что и требовалось.

Для вычисления по этой простой формуле всех <math>n</math> собственных векторов требуется <math>O(n^{2})</math> флопов. К сожалению, формула не обеспечивает численной устойчивости, так как для двух очень близких значений <math>\alpha_{i}</math> может давать неортогональные приближенные собственные векторы <math>u_{i}</math>. Потребовалось целое десятилетие для того, чтобы найти устойчивую альтернативу исходному описанию алгоритма. Снова детали будут обсуждаться позднее в данном разделе.

Алгоритм является рекурсивным.

=== Дефляция ===

До сих пор полагалось, что все <math>d_{i}</math> различны и все <math>u_{i}</math> отличны от нуля. Если это не так, вековое уравнение <math>f(\lambda)=0</math> имеет <math>k</math> вертикальных асимптот, где <math>k<n</math>, а потому <math>k</math> корней. Однако оказывается, что остальные <math>n - k</math> собственных значений могут быть определены без каких-либо усилий: если <math>d_{i}=d_{i+1}</math> или <math>u_{i}=0</math>, то легко показать, что <math>d_{i}</math> является собственным значением и для матрицы <math>D + \rho uu^{T}</math>. В такой ситуации мы говорим о ''дефляции''. На практике выбирается некоторое пороговое значение и дефляция для числа <math>d_{i}</math> регистрируется, если в смысле этого порога <math>d_{i}</math> достаточно близко к <math>d_{i+1}</math> либо <math>u_{i}</math> достаточно мало.

Основной выигрыш от использования дефляции состоит не в том, что убыстряется решение векового уравнения - этот этап в любом случае стоит лишь <math>O(n^{2})</math> операций. Выигрыш заключается в ускорении матричного умножения на последнем шаге алгоритма. Действительно, если <math>u_{i}=0</math>, то соответствующий собственный вектор есть i-й столбец <math>e_{i}</math> единичной матрицы. Это означает, что <math>e_{i}</math> является i-м столбцом в матрице <math>Q_{'}</math>, поэтому при формировании матрицы <math>Q</math> посредством левого умножения <math>Q_{1}</math> на <math>Q_{2}</math> вычисление i-го столбца не требует никаких затрат. Аналогичное упрощение имеет место в случае <math>d_{i} = d_{i+1}</math>. При дефляции многих собственных значений устраняется большая часть работы, связанной с матричным умножением.

== Вычислительное ядро алгоритма ==
Вычислительным ядром последовательной схемы решения является вычисление матрицы <math>Q</math> собственных векторов путём умножения матрицы <math>Q = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}</math> на матрицу <math>Q^{'}</math> Данная операция имеет сложность <math>cn^{3}</math> о чём и говорится в разделе [[#Последовательная сложность алгоритма]] .
Ей предшествует вычисление собственных значений и векторов матрицы <math> D + \rho uu^{T}</math>

== Макроструктура алгоритма ==
В разделе [[#Информационный граф]] описана структура алгоритма, в которой есть блок умножения матриц для вычисления собственных векторов, являющийся вычислительным ядром алгоритма. В соответствующем разделе ([[#Вычислительное ядро алгоритма]]) мы упоминали о том, что данному блоку предшествует вычисление собственных значений, которое производится [https://ru.wikipedia.org/wiki/%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%9D%D1%8C%D1%8E%D1%82%D0%BE%D0%BD%D0%B0 методом Ньютона].

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

Предположим, что некоторое <math>u_{i}</math>, хотя и мало, все же недостаточно мало для того, чтобы была зарегистрирована дефляция. В этом случае применение метода Ньютона к решению векового уравнения встречается с затруднениями. Вспомним, что пересчет приближённого решения <math>u_{j}</math> уравнения <math>f(\lambda) = 0</math> в методе Ньютона основан на следующих положениях:

1. Вблизи точки <math>\lambda = \lambda_{j}</math> функция <math>f(\lambda)</math> аппроксимируется линейной функцией <math>l(\lambda)</math>; график есть прямая линия, касающаяся графика функции <math>f(\lambda)</math> при <math>\lambda = \lambda_{j}</math>.

2. В качестве <math>\lambda_{j+1}</math> берётся нуль этого линейного приближения, т.е. <math>l(\lambda_{j+1})=0</math>.

Функция, показанная на рис.1, не доставляет видимых трудностей методу Ньютона, поскольку вблизи каждого своего нуля <math>f(\lambda)</math> достаточно хорошо аппроксимируется линейными функциями. Однако рассмотрим график функции на рис. 2. Она получена из функции на рис. 1 заменой значения .5 для <math>u_{i}^{2}</math> на .001. Это новое значение недостаточно мало для того, чтобы вызвать дефляцию. График функции в левой части рис.2 визуально не отличим от её вертикальных и горизонтальных асимптот, поэтому в правой части укрупненно воспроизведён фрагмент графика, прилегающий к вертикальной асимптоте <math>\lambda = 2</math>. Видно, что график слишком быстро "выполняет поворот" и для большей части значений <math>\lambda</math> почти горизонтален. Поэтому, применяя метод Ньютона почти к любому начальному приближению <math>\lambda_{0}</math>, мы получаем линейное приближение <math>l(\lambda)</math> с почти горизонтальным графиком и малым положительным угловым коэффициентом. В результате <math>\lambda_{1}</math> является отрицательным числом, огромным по абсолютной величине, которое совершенно бесполезно в качестве приближения к истинному корню.

[[File:GraphNewton.png|thumb|center|800px|Рис. 2. График функции <math> f(\lambda) = 1 + \frac{10^{-3}}{1 - \lambda} + \frac{10^{-3}}{2 - \lambda} + \frac{10^{-3}}{3 - \lambda} + \frac{10^{-3}}{4 - \lambda}</math>]]

Чтобы найти выход из этого положения, можно модифицировать метод Ньютона следующим образом: раз <math>f(\lambda)</math> нельзя хорошо приблизить линейной функцией <math>l(\lambda)</math>, попробуем взять в качестве приближения какую-нибудь другую простую функцию <math>h(\lambda)</math>. Нет ничего особого именно в прямых линиях: для метода Ньютона вместо <math>l(\lambda)</math> можно взять любое приближение <math>h(\lambda)</math>, значения и нули которого легко вычисляются. Функция <math>f(\lambda)</math> имеет полюсы в точках <math>d_{i}</math> и <math>d_{i+1}</math>, которые определяют её поведение в соответствующих окрестностях. Поэтому при поиске корня в интервале <math>(d_{i+1}, d_{i})</math> естественно выбрать функцию <math>h(\lambda)</math>, также имеющую эти полюсы, т.е. функцию вида
<math>h(\lambda)= \frac{c_{1}}{d_{i}-\lambda} + \frac{c_{2}}{d_{i+1}-\lambda} + c_{3}</math>

Константы <math>c_{1},c_{2}</math> и <math>c_{3}</math> обеспечивающие, что <math>h(\lambda)</math> есть приближение к <math>f(\lambda)</math>, можно выбрать несколькими способами. Отметим, что если <math>c_{1},c_{2}</math> и <math>c_{3}</math> уже известны, то уравнение <math>h(\lambda)=0</math> легко решается относительно <math>\lambda</math>, поскольку сводится к эквивалентному квадратному уравнению
<math>c_{1}(d_{i+1}-\lambda)+c_{2}(d_{i}-\lambda)+c_{3}(d_{i}-\lambda)(d_{i+1}-\lambda)=0</math>

Пусть <math>\lambda_{j}</math> - приближённое значение корня. определим <math>c_{1},c_{2}</math> и <math>c_{3}</math> так, чтобы
<math>\frac{c_{1}}{d_{i}-\lambda} + \frac{c_{2}}{d_{i+1}-\lambda} + c_{3} = h(\lambda) \approx f(\lambda) = 1 + \rho \sum_{k=1, n} \frac{u_{k}^{2}} {d_{k}-\lambda} </math>

для <math>\lambda</math> в окрестности <math>\lambda_{j}</math>. Заметим, что

<math>f(\lambda) = 1 + \rho \sum_{k=1, i} \frac{u_{k}^{2}} {d_{k}-\lambda} + \rho \sum_{k=i+1, n} \frac{u_{k}^{2}} {d_{k}-\lambda} \equiv 1 + \psi_{1}(\lambda) + \psi_{2}(\lambda)</math>.

Если <math>\lambda \in (d_{i+1},d_{i})</math>, то <math>\psi_{1}(\lambda)</math> есть сумма положительных слагаемых, а <math>\psi_{2}(\lambda)</math> - сумма отрицательных. Поэтому и <math>\psi_{1}(\lambda)</math>, и <math>\psi_{2}(\lambda)</math> могут быть вычислены с высокой точностью; однако при их сложении вполне вероятно взаимное уничтожение верных разрядов и потеря относительной точности в сумме. Возьмем числа <math>c_{1}</math> и <math>\hat{c_{1}}</math>, такие, что функция

<math>h_{1}(\lambda) \equiv \hat{c_{1}} + \frac{c_{1}}{d_{i}-\lambda}</math> удовлетворяет условиям <math>h_{1}(\lambda_{j}) = \psi_{1}(\lambda_{j})</math> и <math>h_{1}^{'}(\lambda_{j})=\psi_{1}^{'}(\lambda_{j})</math> (1)

Это означает, что гипербола, являющаяся графиком функции <math>h_{1}(\lambda)</math>, касается графика функции <math>\psi_{i}(\lambda)</math> при <math>\lambda = \lambda_{j}</math>. Два условия в (1) - это обычные условия метода Ньютона, за исключением того, что вместо прямой в качестве приближения используется гипербола. Легко проверить, что <math>c_{1}=\psi_{1}^{'}(\lambda_{j})(d_{i} - \lambda_{j})^{2}</math> и <math>\hat{c_{1}}=\psi_{1}(\lambda_{j}) - \psi_{1}^{'}(\lambda_{j})(d_{i} - \lambda_{j})</math>

Подобным же образом выбираем <math>c_{2}</math> и <math>\hat{c_{2}}</math> так, чтобы функция

<math>h_{2}(\lambda) \equiv \hat{c_{2}} + \frac{c_{2}}{d_{i+1}-\lambda}</math>

Удовлетворяла условиям

<math>h_{2}(\lambda_{j}) = \psi_{2}(\lambda_{j})</math> и <math>h_{2}^{'}(\lambda_{j})=\psi_{2}^{'}(\lambda_{j})</math> (2)

Наконец, полагаем

<math>h(\lambda) = 1 + h_{1}(\lambda) + h_{2}(\lambda) = (1 + \hat{c_{1}} + \hat{c_{2}} + \frac{c_{1}}{d_{i}-\lambda} + \frac{c_{2}}{d_{i+1}-\lambda} </math>

<math>\equiv c_{3} + \frac{c_{1}}{d_{i}-\lambda} + \frac{c_{2}}{d_{i+1}-\lambda}</math>.

== Схема реализации последовательного алгоритма ==
''Вычисление собственных значений и собственных векторов симметричной трехдиагональной матрицы посредством стратегии "разделяй и властвуй"'':

''proc dc_eig''<math>(T,Q,\Lambda)...</math> по входной матрице <math>T</math> вычисляются выходные матрицы <math>Q</math> и <math>\Lambda</math>, такие, что <math>T = Q\Lambda Q^{T}</math>

'''Если <math>T</math> - матрица размера <math>1</math> x <math>1</math>'''

1. Присвоить выходным параметрам значения <math> Q = 1, \Lambda = T</math>

'''Иначе'''

1. Представить <math>T</math> в виде <math> T = \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T} </math>

2. ''call dc_eig''<math>(T_{1},Q_{1},\Lambda_{1})</math>

3. ''call dc_eig''<math>(T_{2},Q_{2},\Lambda_{2})</math>

4. Построить <math>D+\rho uu^{T}</math> по <math> \Lambda_{1},\Lambda_{2}, Q_{1}, Q_{2}</math>

5. Найти матрицу собственных значений <math>\Lambda</math>

6. Найти матрицу собственных векторов <math>Q^{'}</math> для матрицы <math>D+\rho uu^{T}</math>

7. Построить матрицу собственных векторов <math>Q</math> для матрицы <math>T</math> : <math> Q = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}* Q^{'} </math>

8. Присвоить выходным параметрам значения <math>Q</math> и <math>\Lambda</math>

'''endif'''

== Последовательная сложность алгоритма ==
Проанализируем сложность алгоритма. Пусть <math>t(n)</math> - число флопов при обработке матрицы размера <math> n x n</math> процедурой ''dc_eig''. Тогда

<math> t(n) = 2t(n/2) </math> два рекурсивных обращения к ''dc_eig''<math>(T_{i},Q_{i},\Lambda_{i})</math>

<math> +O(n^{2})</math> вычисление собственных значений матрицы <math>D+\rho uu^{T}</math>

<math> +O(n^{2})</math> вычисление собственных векторов матрицы <math>D+\rho uu^{T}</math>

<math> +c*n^{3}</math> вычисление матрицы <math>Q = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}*Q^{'}</math>

Если <math> Q_{1}, Q_{2}</math> и <math>Q^{'}</math> рассматриваются как плотные матрицы и используется стандартный алгоритм матричного умножения, то константа <math> c </math> в последней строке равна 1. Таким образом, именно это умножение составляет наиболее трудоёмкую часть алгоритма в целом. Игнорируя члены порядка <math>n^{2}</math>, получаем <math>t(n) = 2t(n/2) + cn^{3}</math>. Решая это разностное уравнение, находим <math> t \approx c\frac{4}{3}n^{3} </math>

На практике константа <math>c</math> обычно гораздо меньше 1, потому что матрица <math>Q^{'}</math> весьма разрежена вследствие явления, называемого дефляцией.

== Информационный граф ==
В данном разделе представлен информационный граф алгоритма: на рисунке 3 изображена часть алгоритма, в которой происходит деление - "Divide", а на рисунке 4 описана часть "Conquer".

[[File:Divide.jpg|thumb|center|left|800px|Рис. 3. Информационный граф алгоритма "Разделяй и властвуй" - шаг Divide]]

[[File:ConquerPart.jpg|thumb|center|1200px|Рис. 4. Информационный граф алгоритма "Разделяй и властвуй" - шаг Conquer]]

Операции:

<math>D</math> - операция Divide разделения матрицы <math> T </math> на матрицы <math>T_{1}</math> и <math>T_{2}</math>

Индексы - размерность

<math>S_{k}</math> - операция Conquer сбора матрицы <math> D </math>

<math> L </math> - операция поиска собственных значений

<math> S^{'} </math> - операция поиска собственных векторов матрицы <math> Q^{'} </math>

<math> S </math> - операция поиска собственных векторов матрицы <math> Q </math>

== Ресурс параллелизма алгоритма ==
Рассмотрев информационный граф алгоритма, можно заметить, что в структуре имеется два параллельных блока - вызов рекурсивных функций ''dc_eig'' для вычисления собственных значений и векторов матриц <math>T_{1}</math> и <math>T_{2}</math> - это единственная часть алгоритма, в которой мы прибегаем к параллелизму. В случае реализации без параллелизма, функции ''dc_eig'' отрабатывают последовательно - сначала со входными параметрами <math>T_{1}, Q_{1}, \Lambda_{1}</math>, затем - <math>T_{2}, Q_{2}, \Lambda_{2}</math>.

== Входные и выходные данные алгоритма ==
На вход алгоритму подаётся трёхдиагональная матрица, описанная в разделе [[#Математическое описание алгоритма]]

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

== Свойства алгоритма ==
Описанный в данной статье алгоритм является самым быстрым алгоритмом, среди существующих: QR / Бисекция и обратная итерация / Разделяй и властвуй.

Последовательная сложность алгоритма: <math> t = c\frac{4}{3}n^{3} </math>

Параллельная сложность алгоритма (в силу использования дефляции): <math> O(n^{2.3)}</math> или в самых редких случаях <math> O(n^{2}) </math> (см. [3] стр. 8).

Особенности алгоритма:

1. Использование дефляции ([[#Дефляция]])

2. Использование адаптированного метода Ньютона ([[#Макроструктура алгоритма]])

= ЧАСТЬ. Программная реализация алгоритма =
Вторая часть описания алгоритмов в рамках AlgoWiki рассматривает все составные части процесса их реализации. Рассматривается как последовательная реализация алгоритма, так и параллельная. Описывается взаимосвязь свойств программ, реализующих алгоритм, и особенностей архитектуры компьютера, на которой они выполняются. Исследуется работа с памятью, локальность данных и вычислений, описывается масштабируемость и эффективность параллельных программ, производительность компьютеров, достигаемая на данной программе. Обсуждаются особенности реализации для разных классов архитектур компьютеров, приводятся ссылки на реализации в существующих библиотеках.

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

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

== Возможные способы и особенности параллельной реализации алгоритма ==
== Масштабируемость алгоритма и его реализации ==
TBD

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

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

== Существующие реализации алгоритма ==
Большого количества реализаций данного алгоритма найдено не было: единственный пример описан на страницах 7, 19 и 20 в [http://www.netlib.org/lapack/lawnspdf/lawn132.pdf этой статье], которая указана в списке литературы.

= Литература =
[1] Дж. Деммель, «Вычислительная линейная алгебра» //С. 230-235

[2] [http://www.netlib.org/utk/people/JackDongarra/PAPERS/104_1999_a-parallel-divide-and-conquer-algorithm.pdf Francoise Tisseury, Jack Dongarra, A Parallel Divide and Conquer Algorithm for the Symmetric Eigenvalue Problem on Distributed Memory Architectures]

[3] [http://www.netlib.org/lapack/lawnspdf/lawn132.pdf Francoise Tisseury, Jack Dongarra, Parallelizing the Divide and Conquer Algorithm for the Symmetric Tridiagonal Eigenvalue Problem on Distributed Memory Architectures]

[4] [https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_algorithm Алгоритм "Разделяй и властвуй" - Wikipedia]

[5] [http://www.cscamm.umd.edu/tadmor/pub/linear-stability/Gill_Tadmor_SISC90.pdf Doron Grill and Eitan Tadmor AN <math>O(N2)</math> METHOD FOR COMPUTING THE EIGENSYSTEM OF <math>N</math>x<math>N</math> SYMMETRIC TRIDIAGONAL MATRICES BY THE DIVIDE AND CONQUER APPROACH]

[[en:Description of algorithm properties and structure]]