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

Основные авторы описания: [[Участник:Роман_Землянский|Р.Землянский]] (Разделы 
[[#Общее описание алгоритма|1.1]], [[#Математическое описание алгоритма|1.2]], [[#Последовательная сложность алгоритма|1.6]], [[#Ресурс параллелизма алгоритма|1.8]], [[#Входные и выходные данные алгоритма|1.9]], [[#Свойства алгоритма|1.10]]), В.Завольсков (Разделы [[#Вычислительное ядро алгоритма|1.3]], [[#Макроструктура алгоритма|1.4]], [[#Схема реализации последовательного алгоритма|1.5]], [[#Информационный граф|1.7]], [[#Существующие реализации алгоритма|2.7]]), [[Участник:IgorS]] (предварительная проверка), [[Участник:Chernyavskiy|А.Ю. Чернявский]] (проверка, редактирование).


{{algorithm
| name              = Метод "Разделяй и властвуй" вычисления собственных значений и векторов трёхдиагональной матрицы
| serial_complexity = <math> c\frac{4}{3}n^{3} </math>
| parallel_complexity = <math> O(n^{2.3}) </math>

крайне редко: <math>  O(n^{2}) </math> 
| pf_width          = 
| input_data        =  <math> 2n + 1 </math>
| output_data       =  <math> n(n + 1) </math>
}}


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

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

:<math>
T = \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} \Lambda_{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}uu^{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} \Lambda_{1} & 0 \\ 0 & \Lambda_{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>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>.

== Вычислительное ядро алгоритма ==
Вычислительным ядром последовательной схемы решения является вычисление матрицы <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 методом Ньютона].

[[File:DivideAndConquerTree.png|thumb|center|left|800px|Рис. 3. Дерево алгоритма "Разделяй и властвуй"]]

[[File:InfoGraph.jpg|thumb|center|1200px|Рис. 4. Детальное описание одного блока алгоритма]]

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


''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> весьма разрежена вследствие явления, называемого дефляцией.

== Информационный граф ==
В данном разделе на рисунке 5 представлен информационный граф алгоритма.

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

Операции: 

<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>.

Таким образом, благодаря параллелизму высота алгоритма параллельного составляет <math>log_{2}(n)</math>, когда как при последовательной реализации высота алгоритма равна сумме геометрической прогрессии <math>1+2+4+8+...+2^{n-1} = 2^{n} - 1</math>, что значительно сокращает число итераций алгоритма.

== Входные и выходные данные алгоритма ==
'''Входные данные:''' 

На вход алгоритму подаётся трёхдиагональная симметричная матрица:

<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}
</math>

'''Объем входных данных:''' 

<math>2n + 1</math> - в силу симметрии, хранить следует только диагональные, поддиагональые и наддиагональные элементы.

'''Выходные данные:'''

На выходе мы получаем вектор собственных значений длины <math> n </math> и  <math> n </math> собственных векторов длины <math> n </math> , соответствующих собственным значениям исходной матрицы.

'''Объем выходных данных:'''

<math> n + n^{2} = n(n + 1) </math>

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

Последовательная сложность алгоритма (не принимая во внимание дефляцию) удовлетворяет рекурсивному соотношению: 

 <math>t(n) = n^3 +2t(\frac{n}{2})</math>

Решение которого имеет вид:

 <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 рассматривает все составные части процесса их реализации. Рассматривается как последовательная реализация алгоритма, так и параллельная. Описывается взаимосвязь свойств программ, реализующих алгоритм, и особенностей архитектуры компьютера, на которой они выполняются. Исследуется работа с памятью, локальность данных и вычислений, описывается масштабируемость и эффективность параллельных программ, производительность компьютеров, достигаемая на данной программе. Обсуждаются особенности реализации для разных классов архитектур компьютеров, приводятся ссылки на реализации в существующих библиотеках.


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

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

== Возможные способы и особенности параллельной реализации алгоритма ==
== Масштабируемость алгоритма и его реализации ==
Исследование проводилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского университета.

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

*	число процессоров [1 : 8] (1,2,3,4,8);
*	размер матрицы [4 : 8192].

Для исследования масштабируемости алгоритм был реализован на библиотеке OpenMP 1.5.5 для языка C++(компилятор g++) на cуперкомпьютере Ломоносов на узле compiler. Для данного узла максимально возможное число потоков - 8. В ходе оценки масштабируемости, программа была запущена в 1, 2, 4, 8 потока. 
Параметры компиляции и запуска: 

 ssh compiler 
 module add slurm
 module add openmpi/1.5.5-icc
 g++ _scratch/test.c++ -fopenmp -std=c++0x -o test 
 sbatch -p test -n1 run ./test

Для управления количеством потоков использовалась операция OpenMP omp_set_num_threads();

[[File:ResultsDiagram.jpg|thumb|center|left|800px|Рис. 6. График времени работы программы в зависимости от числа процессов и размера задачи]]

Масштабируемость рассматриваемой реализации алгоритма '''"Разделяй и властвуй"''':

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

Код исследуемой реализации представлен здесь [http://pastebin.com/QpfYZGkK pastebin.com]

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

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

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

= Литература =
[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]



[[Категория:Статьи в работе]]