Участник:Роман Землянский/Метод "Разделяй и властвуй": различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
 
(не показано 40 промежуточных версий 2 участников)
Строка 1: Строка 1:
Авторы статьи: Завольсков Владислав, Землянский Роман, 614 группа
+
 
 +
Авторы статьи: Завольсков Владислав, [[U:Роман Землянский|Землянский Роман]], 614 группа
  
 
= ЧАСТЬ. Свойства и структура алгоритмов =
 
= ЧАСТЬ. Свойства и структура алгоритмов =
Строка 41: Строка 42:
 
= \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}  
 
= \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}  
  
= \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T}
+
\equiv \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T}
 
</math>
 
</math>
  
Строка 79: Строка 80:
 
Таким образом,
 
Таким образом,
  
<math>det(I + \rho (D - \lambda I)^{-1}uu^{T}) = 1 + \rho u^{T}(D - \lambda I)^{-1}u</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> = 1 + \rho \sum_{i=1, n} \frac{u_{i}^{2}} {d_{i}-\lambda} = 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>).
 
т.е. собственные значения матрицы <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>).
Строка 109: Строка 108:
 
Алгоритм является рекурсивным.
 
Алгоритм является рекурсивным.
  
'''Дефляция'''
+
=== Дефляция ===
  
 
До сих пор полагалось, что все <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>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> достаточно мало.
Строка 116: Строка 115:
  
 
== Вычислительное ядро алгоритма ==
 
== Вычислительное ядро алгоритма ==
Вычислительным ядром последовательной схемы решения является вычисление матрицы <math>Q</math> собственных векторов путём умножения матрицы <math>Q = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}</math> на матрицу <math>Q^{'}</math> Данная операция имеет сложность <math>cn^{3}</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>
 
Ей предшествует вычисление собственных значений и векторов матрицы <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>.
  
 
== Схема реализации последовательного алгоритма ==
 
== Схема реализации последовательного алгоритма ==
Строка 164: Строка 209:
  
  
Если <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 = c\frac{4}{3}n^{3} </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> весьма разрежена вследствие явления, называемого дефляцией.
 
На практике константа <math>c</math> обычно гораздо меньше 1, потому что матрица <math>Q^{'}</math> весьма разрежена вследствие явления, называемого дефляцией.
Строка 170: Строка 215:
 
== Информационный граф ==
 
== Информационный граф ==
 
В данном разделе представлен информационный граф алгоритма: на рисунке 2 изображена структура всего алгоритма в целом, в то время как на рисунке 3 детально описана одна из ячеек структуры.
 
В данном разделе представлен информационный граф алгоритма: на рисунке 2 изображена структура всего алгоритма в целом, в то время как на рисунке 3 детально описана одна из ячеек структуры.
[[File:DivideAndConquerTree.png|thumb|center|left|800px|Рис. 2. Дерево алгоритма "Разделяй и властвуй"]]
+
[[File:DivideAndConquerTree.png|thumb|center|left|800px|Рис. 3. Дерево алгоритма "Разделяй и властвуй"]]
  
[[File:InfoGraph.jpg|thumb|center|1200px|Рис. 3. Детальное описание одного блока алгоритма]]
+
[[File:InfoGraph.jpg|thumb|center|1200px|Рис. 4. Детальное описание одного блока алгоритма]]
  
 
== Ресурс параллелизма алгоритма ==
 
== Ресурс параллелизма алгоритма ==
Строка 178: Строка 223:
  
 
== Входные и выходные данные алгоритма ==
 
== Входные и выходные данные алгоритма ==
На вход алгоритму подаётся трёхдиагональная матрица, описанная в разделе 1.2
+
На вход алгоритму подаётся трёхдиагональная матрица, описанная в разделе [[#Математическое описание алгоритма]]
 +
 
 
На выходе мы получаем собственные значения и собственные вектора исходной матрицы
 
На выходе мы получаем собственные значения и собственные вектора исходной матрицы
  
 
== Свойства алгоритма ==
 
== Свойства алгоритма ==
Какая-нибудь шляпа
+
Описанный в данной статье алгоритм является самым быстрым алгоритмом, среди существующих: 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. Использование адаптированного метода Ньютона ([[#Макроструктура алгоритма]])
  
 
= ЧАСТЬ. Программная реализация алгоритма =
 
= ЧАСТЬ. Программная реализация алгоритма =
Строка 193: Строка 249:
  
 
== Существующие реализации алгоритма ==
 
== Существующие реализации алгоритма ==
TBD
+
Большого количества реализаций данного алгоритма найдено не было: единственный пример описан на страницах 7, 19 и 20 в [http://www.netlib.org/lapack/lawnspdf/lawn132.pdf этой статье], которая указана в списке литературы.
  
 
= Литература =
 
= Литература =
 
[1] Дж. Деммель, «Вычислительная линейная алгебра» //С. 230-235
 
[1] Дж. Деммель, «Вычислительная линейная алгебра» //С. 230-235
  
[2] Francoise Tisseury, Jack Dongarra, Parallelizing the Divide and Conquer Algorithm for the Symmetric Tridiagonal Eigenvalue Problem on Distributed Memory Architectures
+
[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] Francoise Tisseury, Jack Dongarra, A Parallel Divide and Conquer Algorithm for the Symmetric Eigenvalue Problem on Distributed Memory Architectures
 
  
[4] Per-Olof Persson, Lecture 17 Other Eigenvalue Algorithms //
+
[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]
  
[5] [https://www.cscamm.umd.edu%2Ftadmor%2Fpub%2Flinear-stability%2FGill_Tadmor_SISC90.pdf]
+
[4] [https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_algorithm Алгоритм "Разделяй и властвуй" - Wikipedia]
  
[6][https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_algorithm]
+
[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]]
 
[[en:Description of algorithm properties and structure]]

Текущая версия на 23:49, 1 ноября 2016

Авторы статьи: Завольсков Владислав, Землянский Роман, 614 группа

1 ЧАСТЬ. Свойства и структура алгоритмов

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

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

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

Пусть

[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} \lt = \ldots \lt =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} \lt \gt 0[/math] (случай общего положения), то [math]f(\lambda)[/math] имеет график типа показанного на рис.1(где [math]n = 4[/math] и [math]\rho \gt 0[/math]).

Рис. 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}}\gt 0 [/math], функция возрастает всюду, кроме точек [math]\lambda = d_{i}[/math]. Поэтому корни функции разделяются числами [math]d_{i}[/math] и ещё один корень находится справа от точки [math]d_{1}[/math] (на рис. 1 [math]d_{1} = 4[/math]). (При [math]\rho\lt 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]. Потребовалось целое десятилетие для того, чтобы найти устойчивую альтернативу исходному описанию алгоритма. Снова детали будут обсуждаться позднее в данном разделе.

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

1.2.1 Дефляция

До сих пор полагалось, что все [math]d_{i}[/math] различны и все [math]u_{i}[/math] отличны от нуля. Если это не так, вековое уравнение [math]f(\lambda)=0[/math] имеет [math]k[/math] вертикальных асимптот, где [math]k\lt 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]. При дефляции многих собственных значений устраняется большая часть работы, связанной с матричным умножением.

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

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

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

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

1.4.1 Решение векового уравнения

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

Предположим, что некоторое [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] является отрицательным числом, огромным по абсолютной величине, которое совершенно бесполезно в качестве приближения к истинному корню.

Рис. 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].

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

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


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

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

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

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

В данном разделе представлен информационный граф алгоритма: на рисунке 2 изображена структура всего алгоритма в целом, в то время как на рисунке 3 детально описана одна из ячеек структуры.

Рис. 3. Дерево алгоритма "Разделяй и властвуй"
Рис. 4. Детальное описание одного блока алгоритма

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

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

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

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

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

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

Описанный в данной статье алгоритм является самым быстрым алгоритмом, среди существующих: 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. Использование адаптированного метода Ньютона (#Макроструктура алгоритма)

2 ЧАСТЬ. Программная реализация алгоритма

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


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

TBD


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

Большого количества реализаций данного алгоритма найдено не было: единственный пример описан на страницах 7, 19 и 20 в этой статье, которая указана в списке литературы.

3 Литература

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

[2] Francoise Tisseury, Jack Dongarra, A Parallel Divide and Conquer Algorithm for the Symmetric Eigenvalue Problem on Distributed Memory Architectures

[3] Francoise Tisseury, Jack Dongarra, Parallelizing the Divide and Conquer Algorithm for the Symmetric Tridiagonal Eigenvalue Problem on Distributed Memory Architectures

[4] Алгоритм "Разделяй и властвуй" - Wikipedia

[5] 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