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

Материал из Алговики
Перейти к навигации Перейти к поиску
(Новая страница: «Метод «разделяй и властвуй» вычисления собственных значений и векторов симметричной тр…»)
 
 
(не показано 38 промежуточных версий 4 участников)
Строка 1: Строка 1:
 +
{{Assignment|ASA}}
 +
 
Метод «разделяй и властвуй» вычисления собственных значений и векторов симметричной трехдиагональной матрицы.
 
Метод «разделяй и властвуй» вычисления собственных значений и векторов симметричной трехдиагональной матрицы.
−
Автор: Изотова И.А. (619 группа)
+
 
 +
Авторы:
 +
 +
Изотова И.А. (619 группа) пункты 1.1, 1.2, 1.5, 1.7, 1.9
 +
 
 +
Ульянов Г.В. (619 группа) пункты 1.3, 1.4, 1.6, 1.8, 1.10
  
  
Строка 6: Строка 13:
  
 
== Общее описание алгоритма ==
 
== Общее описание алгоритма ==
−
Метод «разделяй и властвуй» - класс алгоритмов вычисления собственных значений и векторов симметричной трехдиагональной матрицы. Наиболее быстрый из существующих методов вычисления собственных значений и векторов симметричной трехдиагональной матрицы начиная с порядка <math>n</math>, примерно равного 26. (Точное значение этого порогового порядка зависит от компьютера.) Его численно устойчивая реализация весьма не тривиальна. В самом деле, хотя впервые метод был предложен еще в 1981 г., "правильный" способ его реализации был найден лишь в 1992 г. Этот способ воплощен LAPACK-программами ssyevd (для плотных матриц) и sstevd (для трехдиагональных матриц). В них стратегия "разделяй-и-влавствуй" используется для матриц порядка, большего чем 25. Для матриц меньшего порядка (или если нужны только собственные значения) происходит автоматический переход к QR-итерации.  
+
'''«Разделяй и властвуй»''' - метод вычисления собственных значений и векторов симметричной трехдиагональной матрицы. В настоящее время это один из самых быстрых методов при вычислении собственных значений и векторов симметричной трехдиагональной матрицы порядка <math>n > 25</math> (точное значение этого порогового порядка зависит от компьютера.) Однако алгоритм не гарантирует высокой относительной точности для малых собственных значений.Его численно устойчивая реализация весьма не тривиальна. Хотя впервые метод был предложен еще в 1981 г., "правильный" способ его реализации был найден лишь в 1992 г. Этот способ воплощен LAPACK-программами ssyevd (для плотных матриц) и sstevd (для трехдиагональных матриц). В них стратегия "разделяй-и-влавствуй" используется для матриц порядка, большего чем 25. Для матриц меньшего порядка (или если нужны только собственные значения) происходит автоматический переход к QR-итерации. Однако, из-за высокой сложности, алгоритм пока не вошел в состав ALGLIB.
 +
 
 +
В наихудшем случае алгоритм «разделяй-и-властвуй» требует <math> O(n^{3})</math> операций, но константа, на практике оказывается весьма малой. Для некоторых специальных распределений собственных значений — <math> O(n^{2})</math> (подробнее в разделе [[#Последовательная сложность алгоритма|1.6 Последовательная сложность алгоритма]]).
  
 
== Математическое описание алгоритма ==
 
== Математическое описание алгоритма ==
−
==Разделяй==
+
'''Разделяй'''
  
−
Пусть матрицы <math>T</math>, <math>T_{1}</math>, <math>T_{2}</math> имеют размерность <math>n \times n</math>, <math>m \times m</math>, <math>(n - m) \times (n - m)</math>.
+
Пусть матрицы <math>T, T_{1}</math>, <math>T_{2}</math> имеют размерность <math>n \times n</math>, <math>m \times m</math>, <math>(n - m) \times (n - m)</math>.
  
 
Запишем матрицу <math>T</math> в следующем виде:
 
Запишем матрицу <math>T</math> в следующем виде:
Строка 53: Строка 62:
  
 
Предположим, что нам известны спектральные разложения матриц <math>\hat{T}_{1}</math> и <math>\hat{T}_{2}</math>, т.е.:
 
Предположим, что нам известны спектральные разложения матриц <math>\hat{T}_{1}</math> и <math>\hat{T}_{2}</math>, т.е.:
−
[[diagonalizable matrix|diagonalization]]s <math>\hat{T}_{1} = Q_{1} D_{1} Q_{1}^{T}</math> и <math>\hat{T}_{2} = Q_{2} D_{2} Q_{2}^{T}</math>.   
+
<math>\hat{T}_{1} = Q_{1} D_{1} Q_{1}^{T}</math> и <math>\hat{T}_{2} = Q_{2} D_{2} Q_{2}^{T}</math>.   
  
−
==Властвуй==
+
'''Властвуй'''
  
 
Установим связь между собственными значениями матрицы <math>T</math> и собственными значениями матриц <math>T_{1}</math> и <math>T_{2}</math>:
 
Установим связь между собственными значениями матрицы <math>T</math> и собственными значениями матриц <math>T_{1}</math> и <math>T_{2}</math>:
Строка 62: Строка 71:
 
T = \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T}  
 
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} \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}
+
= \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>,
  
 
где <math> u = \begin{bmatrix} Q_{1}^{T} & 0 \\ 0 & Q_{2}^{T}\end{bmatrix}v </math> - матрица, в первой строке которой находятся элементы последнего столбца матрицы <math> Q_{1}^{T}</math>, а во второй - элементы первого столбца матрицы <math> Q_{2}^{T}</math>, так как <math>v = \begin{bmatrix} 0 , \ldots  ,  0 , 1 , 1 , 0 \ldots , 0 \end{bmatrix}^T</math>.
 
где <math> u = \begin{bmatrix} Q_{1}^{T} & 0 \\ 0 & Q_{2}^{T}\end{bmatrix}v </math> - матрица, в первой строке которой находятся элементы последнего столбца матрицы <math> Q_{1}^{T}</math>, а во второй - элементы первого столбца матрицы <math> Q_{2}^{T}</math>, так как <math>v = \begin{bmatrix} 0 , \ldots  ,  0 , 1 , 1 , 0 \ldots , 0 \end{bmatrix}^T</math>.
 
Следовательно, <math>T</math> имеет такие же собственные значения, что и подобная ей матрица <math>D + \rho uu^{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>D =\begin{bmatrix} \Lambda_{1} & 0 \\ 0 & \Lambda_{2}\end{bmatrix}</math> - диагональная матрица,  
 
<math>\rho = b_{m}</math> - число, а <math>u</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_{1}, \ldots, d_{n}</math> матрицы <math>D</math> упорядочены по не возрастанию: <math>d_{n} <= \ldots <=d_{1}</math>.
Строка 84: Строка 93:
 
т.е. собственные значения матрицы <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>).
  
−
[[File:График_функции.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>]]
+
[[File:График_функции.PNG|thumb|center|600px|Рис. 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>]]
  
 
По Рис.1 можно заметить, что прямая <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> флопов.
 
По Рис.1 можно заметить, что прямая <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> флопов.
Строка 91: Строка 100:
 
   ''Лемма 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> флопов.
 
   ''Лемма 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] =  
 
<math>(D + \rho uu^{T})[(D - \alpha I)^{-1}u] =  
  
Строка 106: Строка 116:
 
Что и требовалось доказать.
 
Что и требовалось доказать.
  
−
Получается, что для вычисления по этой простой формуле всех <math>n</math> собственных векторов потребуется <math>O(n^{2})</math> флопов. К сожалению, данная формула не обеспечивает нам численной устойчивости, так как для двух очень близких значений <math>\alpha_{i}</math> может давать неортогональные приближенные собственные векторы <math>u_{i}</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>Q</math> собственных векторов при помощи умножения матрицы <math>Q = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}</math> на матрицу <math>Q^{'}</math>. Эта операция имеет сложность <math>cn^{3}</math> о чём будет говориться ниже, в разделе [[#Последовательная сложность алгоритма|1.6 Последовательная сложность алгоритма]].
−
Вычислительным ядром последовательной схемы решения будет являться вычисление матрицы <math>Q</math> собственных векторов при помощи умножения матрицы <math>Q = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}</math> на матрицу <math>Q^{'}</math>. Эта операция имеет сложность <math>cn^{3}</math> о чём будет говориться ниже, в разделе [[1.6 Последовательная сложность алгоритма]] .
 
  
 
Данной операции предшествует вычисление собственных значений и векторов матрицы <math> D + \rho uu^{T}</math>.
 
Данной операции предшествует вычисление собственных значений и векторов матрицы <math> D + \rho uu^{T}</math>.
Строка 116: Строка 129:
 
== Макроструктура алгоритма ==
 
== Макроструктура алгоритма ==
  
−
В разделе [[#Информационный граф]] описана структура алгоритма, в которой есть блок умножения матриц для вычисления собственных векторов, являющийся вычислительным ядром алгоритма. В соответствующем разделе ([[#Вычислительное ядро алгоритма]]) мы упоминали о том, что данному блоку предшествует вычисление собственных значений, которое производится [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 методом Ньютона].
+
Алгоритм "разделяй и властвуй" для вычисления собственных значений и векторов симметричной трехдиагональной матрицы состоит из вызова рекурсивной функции ''proc dc_eig'' <math>(T, Q, \Lambda)</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:Рисунок2.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>T</math> представляется в виде <math> T = \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T} </math>, представление в таком виде возможно  до тех пор, пока размерность матрицы <math>T</math> не достгнет 1.  
−
<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>T_{1}</math> и <math>T_{2}</math> также вызываются соотвтетствующие функции ''proc dc_eig'' <math>(T_{1}, Q_{1}, \Lambda_{1})</math> и ''proc dc_eig'' <math>(T_{2}, Q_{2}, \Lambda_{2})</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> (*)
 
  
−
Это означает, что гипербола, являющаяся графиком функции <math>h_{1}(\lambda)</math>, касается графика функции <math>\psi_{i}(\lambda)</math> при <math>\lambda = \lambda_{j}</math>. Два условия в (*) - это обычные условия метода Ньютона, за исключением того, что вместо прямой в качестве приближения используется гипербола. Легко проверить, что <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>D+\rho uu^{T}</math> по <math> \Lambda_{1},\Lambda_{2}, Q_{1}, Q_{2}</math>.
  
−
Подобным же образом выбираем <math>c_{2}</math> и <math>\hat{c_{2}}</math> так, чтобы функция
+
*Вычисление собственных значений<math>\Lambda</math> и векторов <math>Q^{'}</math> для матрицы <math>D+\rho uu^{T}</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>
 
  
−
Наконец, полагаем
+
*Построение матрицы собственных векторов <math>Q</math>для матрицы <math>T</math>.
−
<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} \equiv c_{3} + \frac{c_{1}}{d_{i}-\lambda} + \frac{c_{2}}{d_{i+1}-\lambda}</math>.
 
−
 
 
−
Как только из векового уравнения найдены собственные значения alpha I матрицы <math>D + \rho uu^{T}</math>, можно вычислить собственные векторы, пользуясь просто формулой <math>(D - \alpha I)^{-1}u</math> и из Леммы 2.К сожалению, вычисления по этой формуле могут быть неустойчивы. Так будет, в частности, когда собственные значения alpha I, alpha_{i+1} очень близки. Интуитивно проблема состоит в том, что выражения <math>(D - \alpha I)^{-1}u</math> и <math>(D - \alpha_{i+1})^{-1}u</math> "очень похожи", в то время как требуется получить ортогональные собственные векторы. Более точно, если alpha I и alpha_{i+1} очень близки, оба близки к находящемуся между ними числу d_{i}. Поэтому имеет место заметная потеря верных знаков при вычислении d_{i}-alpha I или alpha_{i+1}-d_{i}, либо при решении векового уравнения методом Ньютона. В любом из этих случаев d_{i}-alpha_{i} и alpha_{i+1}1-d_{i} могут содержать большие относительные ошибки, вследствие чего вычисленные собственные векторы <math>(D - \alpha I)^{-1}u</math> и <math>(D - \alpha_{i+1})^{-1}u</math> и очень неточны, и далеки от ортогональности.
 
  
 
== Схема реализации последовательного алгоритма ==
 
== Схема реализации последовательного алгоритма ==
−
Рассмотри реализацию метода "разделяй и властвуй" в виде алгоритма:
+
Рассмотри реализацию метода "разделяй и властвуй" в виде алгоритма, состоящего из частного случая, когда матрица <math>T</math> имеет размерность <math>1 \times 1</math> и общего случая, когда матрица <math>T</math> - симметричная трехдиагональная матрица произвольной размерности, кроме <math>1 \times 1</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>
+
  ''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 \times 1</math>
+
     ''if'' матрица <math>T</math> имеет размерность <math>1 \times 1</math>                       (''проверяем матрицу <math>T</math> на принадлежность частному случаю алгоритма'')
−
         присвоить выходным параметрам значения <math>Q = 1, \Lambda = T</math>
+
    ''then''
−
     иначе
+
         присвоить выходным параметрам значения <math>Q = 1, \Lambda = T</math>             (''если матрица <math>T</math> <math>1 \times 1</math>, то СЗ и СВ матрицы <math>T</math> найдены'')
−
         <math>T</math> в виде <math> T = \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T} </math>
+
     ''else''                                                                                    (''если матрица <math>T</math> имеет общей вид, то:'')
−
         ''call dc_eig''<math>(T_{1},Q_{1},\Lambda_{1})</math>
+
         <math>T</math> в виде <math> T = \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T} </math>                                                                 (''разложим матрицу <math>T</math>'')
−
         ''call dc_eig''<math>(T_{2},Q_{2},\Lambda_{2})</math>
+
         ''call dc_eig''<math>(T_{1},Q_{1},\Lambda_{1})</math>                             (''вызываем рекурсивно функцию ''dc_eig'', вычисляющую СЗ и СВ матрицы <math>T_1</math>'')
−
         <math>D+\rho uu^{T}</math> по <math> \Lambda_{1},\Lambda_{2}, Q_{1}, Q_{2}</math>
+
         ''call dc_eig''<math>(T_{2},Q_{2},\Lambda_{2})</math>                             (''вызываем рекурсивно функцию ''dc_eig'', вычисляющую СЗ и СВ матрицы <math>T_2</math>'')
 +
         Построить <math>D+\rho uu^{T}</math> по <math> \Lambda_{1},\Lambda_{2}, Q_{1}, Q_{2}</math>
 
         Найти матрицу собственных значений <math>\Lambda</math>
 
         Найти матрицу собственных значений <math>\Lambda</math>
 
         Найти матрицу собственных векторов <math>Q^{'}</math> для матрицы <math>D+\rho uu^{T}</math>
 
         Найти матрицу собственных векторов <math>Q^{'}</math> для матрицы <math>D+\rho uu^{T}</math>
Строка 174: Строка 161:
  
 
== Последовательная сложность алгоритма ==
 
== Последовательная сложность алгоритма ==
−
В данном разделе описания свойств алгоритма приводится оценка его [[глоссарий#Последовательная сложность|''последовательной сложности'']], т.е. числа операций, которые нужно выполнить при последовательном исполнении алгоритма (в соответствии с [[#Описание схемы реализации последовательного алгоритма|п.1.5]]). Для разных алгоритмов понятие операции, в терминах которой оценивается его сложность, может существенно различаться. Это могут быть операции для работы с вещественными числами, целыми числами, поразрядные операции, обращения в память, обновления элементов массива, элементарные функции, макрооперации и другие. В LU-разложении преобладают арифметические операции над вещественными числами, а для транспонирования матриц важны лишь обращения к памяти: это и должно найти отражение в описании.
+
Пусть <math>t(n)</math> - число флопов при обработке матрицы размера <math>n \times n</math> процедурой ''dc_eig''. Тогда
  
−
Если выбор конкретного типа операций для оценки сложности алгоритма не очевиден, то нужно привести обоснование возможных вариантов. В некоторых случаях можно приводить оценку не всего алгоритма, а лишь его вычислительного ядра: в таком случае это нужно отметить, сославшись [[#Общее описание алгоритма|на п.1.1]].
+
<math>t(n) = 2t(n/2)</math> два рекурсивных обращения к ''dc_eig''<math>(T_{i},Q_{i},\Lambda_{i})</math>
  
−
Например, сложность алгоритма суммирования элементов вектора сдваиванием равна <math>n-1</math>. Сложность быстрого преобразования Фурье (базовый алгоритм Кули-Тьюки) для векторов с длиной, равной степени двойки – <math>n\log_2n</math> операций комплексного сложения и <math>(n\log_2n)/2</math> операций комплексного умножения. Сложность базового алгоритма разложения Холецкого (точечный вариант для плотной симметричной и положительно-определенной матрицы) это <math>n</math> вычислений квадратного корня, <math>n(n-1)/2</math> операций деления, по <math>(n^3-n)/6</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<<1</math> ( так как матрица <math>Q^{'}</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 = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}*Q^{'}</math> является самым затратным действием, то дефляция многих собственных значений может привести в снижению сложности до <math>O(n^{2})</math> операций.
  
 
== Информационный граф ==
 
== Информационный граф ==
−
Это очень важный раздел описания. Именно здесь можно показать (увидеть) как устроена параллельная структура алгоритма, для чего приводится описание и изображение его информационного графа ([[глоссарий#Граф алгоритма|''графа алгоритма'']] [1]). Для рисунков с изображением графа будут составлены рекомендации по их формированию, чтобы все информационные графы, внесенные в энциклопедию, можно было бы воспринимать и интерпретировать одинаково. Дополнительно можно привести полное параметрическое  описание графа в терминах покрывающих функций [1].
+
В данном разделе представлен информационный граф алгоритма: на рисунке 2 изображена часть алгоритма, в которой происходит деление - "разделяй", а на рисунке 3 описана часть "властвуй".
 +
 
 +
[[File:Divide.png|thumb|center|left|8000px|Рис. 2 Информационный граф шага "разделяй"]]
  
−
Интересных вариантов для отражения информационной структуры алгоритмов много. Для каких-то алгоритмов нужно показать максимально подробную структуру, а иногда важнее макроструктура. Много информации несут разного рода проекции информационного графа, выделяя его регулярные составляющие и одновременно скрывая несущественные детали. Иногда оказывается полезным показать последовательность в изменении графа при изменении значений внешних переменных  (например, размеров матриц): мы часто ожидаем "подобное" изменение информационного графа, но это изменение не всегда очевидно на практике. 
+
[[File:Divide and Conq.png|thumb|center|left|8000px|Рис. 3 Связь шагов "Разделяй" и "властвуй"]]
  
−
В целом, задача изображения графа алгоритма весьма нетривиальна. Начнем с того, что это потенциально бесконечный граф, число вершин и дуг которого определяется значениями внешних переменных, а они могут быть весьма и весьма велики. В такой ситуации, как правило, спасают упомянутые выше соображения подобия, делающие графы для разных значений внешних переменных "похожими": почти всегда достаточно привести лишь один граф небольшого размера, добавив, что графы для остальных значений будут устроены "точно также". На практике, увы, не всегда все так просто, и здесь нужно быть аккуратным.
+
[[File:Conq.png|thumb|center|left|8000px|Рис. 4 Информационный граф шага "властвуй"]]
  
−
Далее, граф алгоритма - это потенциально многомерный объект. Наиболее естественная система координат для размещения вершин и дуг информационного графа опирается на структуру вложенности циклов в реализации алгоритма. Если глубина вложенности циклов не превышает трех, то и граф размещается в привычном трехмерном пространстве, однако для более сложных циклических конструкций с глубиной вложенности 4 и больше необходимы специальные методы представления и изображения графов.
+
Операции:
  
−
В данном разделе AlgoWiki могут использоваться многие интересные возможности, которые еще подлежат обсуждению: возможность повернуть граф при его отображении на экране компьютера для выбора наиболее удобного угла обзора, разметка вершин по типу соответствующим им операций, отражение [[глоссарий#Ярусно-параллельная форма графа алгоритма|''ярусно-параллельной формы графа'']] и другие. Но в любом случае нужно не забывать главную задачу данного раздела - показать информационную структуру алгоритма так, чтобы стали понятны все его ключевые особенности, особенности параллельной структуры, особенности множеств дуг, участки регулярности и, напротив, участки с недерминированной структурой, зависящей от входных данных.
+
*<math>D</math> - операция "разделяй" разделения матрицы <math> T </math> на матрицы <math>T_{1}</math> и <math>T_{2}</math>
  
−
На рис.1 показана информационная структура алгоритма умножения матриц, на рис.2 - информационная структура одного из вариантов алгоритма решения систем линейных алгебраических уравнений с блочно-двухдиагональной матрицей.
+
*<math>C_{n}</math> - операция "властвуй" сбора матрицы <math> D_{n*n} </math>
  
−
[[file:Fig1.svg|thumb|center|300px|Рис.1. Информационная структура алгоритма умножения матриц]]
+
*<math> R_{n} </math> - операция поиска собственных значений матрицы <math> T_{n*n} </math>
−
[[file:Fig2.svg|thumb|center|300px|Рис.2. Информационная структура одного из вариантов алгоритма решения систем линейных алгебраических уравнений с блочно-двухдиагональной матрицей]]
 
  
−
== Ресурс параллелизма алгоритма ==
+
*<math> S^{'}_{n} </math> - операция поиска собственных векторов матрицы <math> D_{n*n} </math>
−
Здесь приводится оценка [[глоссарий#Параллельная сложность|''параллельной сложности'']] алгоритма: числа шагов, за которое можно выполнить данный алгоритм в предположении доступности неограниченного числа необходимых процессоров (функциональных устройств, вычислительных узлов, ядер и т.п.). Параллельная сложность алгоритма понимается как высота канонической ярусно-параллельной формы [1]. Необходимо указать, в терминах каких операций дается оценка. Необходимо описать сбалансированность параллельных шагов по числу и типу операций, что определяется шириной ярусов канонической ярусно-параллельной формы и составом операций на ярусах.
 
  
−
Параллелизм в алгоритме часто имеет естественную иерархическую структуру. Этот факт очень полезен на практике, и его необходимо отразить в описании. Как правило, подобная иерархическая структура параллелизма хорошо отражается в последовательной реализации алгоритма через циклический профиль результирующей программы (конечно же, с учетом графа вызовов), поэтому циклический профиль ([[#Описание схемы реализации последовательного алгоритма|п.1.5]]) вполне  может быть использован и для отражения ресурса параллелизма.
+
*<math> S_{n} </math> - операция поиска собственных векторов матрицы <math> T_{n*n} </math>
  
−
Для описания ресурса параллелизма алгоритма (ресурса параллелизма информационного графа) необходимо указать ключевые параллельные ветви в терминах [[глоссарий#Конечный параллелизм|''конечного'']] и [[глоссарий#Массовый параллелизм|''массового'']] параллелизма. Далеко не всегда ресурс параллелизма выражается просто, например, через [[глоссарий#Кооодинатный параллелизм|''координатный параллелизм'']] или, что то же самое, через независимость итераций некоторых циклов (да-да-да, циклы - это понятие, возникающее лишь на этапе реализации, но здесь все так связано… В данном случае, координатный параллелизм означает, что информационно независимые вершины лежат на гиперплоскостях, перпендикулярных одной из координатных осей). С этой точки зрения, не менее важен и ресурс [[глоссарий#Скошенный параллелизм|''скошенного параллелизма'']]. В отличие от координатного параллелизма, скошенный параллелизм намного сложнее использовать на практике, но знать о нем необходимо, поскольку иногда других вариантов и не остается: нужно оценить потенциал алгоритма, и лишь после этого, взвесив все альтернативы, принимать решение о конкретной параллельной реализации. Хорошей иллюстрацией может служить алгоритм, структура которого показана на рис.2: координатного параллелизма нет, но есть параллелизм скошенный, использование которого снижает сложность алгоритма с <math>n\times m</math> в последовательном случае до <math>(n+m-1)</math> в параллельном варианте.
+
Данные:
  
−
Рассмотрим алгоритмы, последовательная сложность которых уже оценивалась в [[#Последовательная сложность алгоритма|п.1.6]]. Параллельная сложность алгоритма суммирования элементов вектора сдваиванием равна <math>\log_2n</math>, причем число операций на каждом ярусе убывает с <math>n/2</math> до <math>1</math>. Параллельная сложность быстрого преобразования Фурье (базовый алгоритм Кули-Тьюки) для векторов с длиной, равной степени двойки - <math>\log_2n</math>. Параллельная сложность базового алгоритма разложения Холецкого (точечный вариант для плотной симметричной и положительно-определенной матрицы) это <math>n</math> шагов для вычислений квадратного корня, <math>(n-1)</math> шагов для операций деления и <math>(n-1)</math> шагов для операций умножения и сложения.
+
*<math>Q_{n}</math> - вектор собственных значений размерности <b>n</b>.  
  
−
== Входные и выходные данные алгоритма ==
+
*<math>V_{n*n}</math> - матрица собственных векторов размерности <b>n</b>.
−
В данном разделе необходимо описать объем, структуру, особенности и свойства входных и выходных данных алгоритма: векторы, матрицы, скаляры, множества, плотные или разреженные структуры данных, их объем. Полезны предположения относительно диапазона значений или структуры, например, диагональное преобладание в структуре входных матриц, соотношение между размером матриц по отдельным размерностям, большое число матриц очень малой размерности, близость каких-то значений к машинному нулю, характер разреженности матриц и другие.
 
  
−
== Свойства алгоритма ==
+
== Ресурс параллелизма алгоритма ==
−
Описываются прочие свойства алгоритма, на которые имеет смысл обратить внимание на этапе реализации. Как и ранее, никакой привязки к конкретной программно-аппаратной платформе не предполагается, однако вопросы реализации в проекте AlgoWiki всегда превалируют, и необходимость обсуждения каких-либо свойств алгоритмов определяется именно этим.
+
Параллелизм достигается за счёт вызова рекурсивных функций ''dc_eig'', вычисляющих собственные значения и собственные вектора матриц <math>T_1</math> и <math>T_2</math>. Рекурсивное разделение матрицы приводит к возможности построения структуры в виде двоичного дерева, в каждой вершине которого происходит вычисение функции ''dc_eig''. Таким образом высота ярусно-параллельной формы, где вычисления для каждого листа дерева в такой структуре происходят независимо, составляет <math>log_2(n)</math>.
  
−
Весьма полезным является ''соотношение последовательной и параллельной сложности'' алгоритма. Оба понятия мы рассматривали ранее, но здесь делается акцент на том выигрыше, который теоретически может дать параллельная реализация алгоритма. Не менее важно описать и те сложности, которые могут возникнуть в процессе получения параллельной версии алгоритма.
+
== Входные и выходные данные алгоритма ==
 +
'''Входные данные:'''
  
−
[[глоссарий#Вычислительная мощность|''Вычислительная мощность'']] алгоритма равна отношению числа операций к суммарному объему входных и выходных данных. Она показывает, сколько операций приходится на единицу переданных данных. Несмотря на простоту данного понятия, это значение исключительно полезно на практике: чем выше  вычислительная мощность, тем меньше накладных расходов вызывает перемещение данных для их обработки, например, на сопроцессоре, ускорителе или другом узле кластера. Например, вычислительная мощность скалярного произведения двух векторов равна всего лишь <math>1</math>, а вычислительная мощность алгоритма умножения двух квадратных матриц равна <math>2n/3</math>.
+
Матрица <math>T</math> размерностью <math>n \times n</math>:
 +
:<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}
 +
</math>
  
−
Вопрос первостепенной важности на последующем этапе реализации - это [[глоссарий#Устойчивость|''устойчивость'']] алгоритма. Все, что касается различных сторон этого понятия, в частности, оценки устойчивости, должно быть описано в данном разделе.
+
'''Объем входных данных:'''
  
−
''Сбалансированность'' вычислительного процесса можно рассматривать с разных сторон. Здесь и сбалансированность типов операций, в частности, арифметических операций между собой (сложение, умножение, деление) или же арифметических операций по отношению к операциям обращения к памяти (чтение/запись). Здесь и сбалансированность операций между параллельными ветвями алгоритма. С одной стороны, балансировка нагрузки является необходимым условием эффективной реализации алгоритма. Вместе с этим, это очень непростая задача, и в описании должно быть отмечено явно, насколько алгоритм обладает этой особенностью. Если обеспечение сбалансированности не очевидно, желательно описать возможные пути решения этой задачи.
+
<math>2n-1</math> - хранить следует только диагональные, поддиагональные и наддиагональные элементы.
  
−
На практике важна [[глоссарий#Детерминированность|''детерминированность алгоритмов'']], под которой будем понимать постоянство структуры вычислительного процесса. С этой точки зрения, классическое умножение плотных матриц является детерминированным алгоритмом, поскольку его структура при фиксированном размере матриц никак не зависит от элементов входных матриц. Умножение разреженных матриц, когда матрица хранятся в одном из специальных форматов, свойством детерминированности уже не обладает: его свойства, например, степень локальности данных зависит от структуры разреженности входных матриц. Итерационный алгоритм с выходом по точности также не является детерминированным: число итераций, а значит и число операций, меняется в зависимости от входных данных. В этом же ряду стоит использование датчиков случайных чисел, меняющих вычислительный процесс для различных запусков программы. Причина выделения свойства детерминированности понятна: работать с детерминированным алгоритмом проще, поскольку один раз найденная структура и будет определять качество его реализации. Если детерминированность нарушается, то это должно быть здесь описано вместе с описанием того, как недетерминированность влияет на структуру вычислительного процесса.
+
'''Выходные данные:'''
  
−
Серьезной причиной недетерминированности работы параллельных программ является изменение порядка выполнения ассоциативных операций. Типичный пример - это использование глобальных MPI-операций на множестве параллельных процессов, например, суммирование элементов распределенного массива. Система времени исполнения MPI сама выбирает порядок выполнения операций, предполагая выполнение свойства ассоциативности, из-за чего ошибки округления меняются от запуска программы к запуску, внося изменения в конечный результат ее работы. Это очень серьезная проблема, которая сегодня встречается часто на системах с массовым параллелизмом и определяет отсутствие повторяемости результатов работы параллельных программ. Данная особенность характерна для [[#ЧАСТЬ. Программная реализация алгоритмов|второй части AlgoWiki]], посвященной реализации алгоритмов, но вопрос очень важный, и соответствующие соображения, по возможности, должны быть отмечены и здесь.
+
Собственные значения <math>\alpha_{i}</math> (всего <math>n</math>) и собственные векторы <math>\lambda_{i}</math> (всего <math>n^2</math>).
  
−
Заметим, что, в некоторых случаях, недетерминированность в структуре алгоритмов можно "убрать" введением соответствующих макроопераций, после чего структура становится не только детерминированной, но и более понятной для восприятия. Подобное действие также следует отразить в данном разделе.
+
'''Объем выходных данных:'''
  
−
[[глоссарий#Степень исхода|''Степень исхода вершины информационного графа'']] показывает, в скольких операциях ее результат будет использоваться в качестве аргумента. Если степень исхода вершины велика, то на этапе реализации алгоритма нужно позаботиться об эффективном доступе к результату ее работы. В этом смысле, особый интерес представляют рассылки данных, когда результат выполнения одной операции используется во многих других вершинах графа, причем число таких вершин растет с увеличением значения внешних переменных.
+
<math>n+n^2</math>
  
−
''"Длинные" дуги в информационном графе'' [1] говорят о потенциальных сложностях с размещением данных в иерархии памяти компьютера на этапе выполнения программы. С одной стороны, длина дуги зависит от выбора конкретной системы координат, в которой расположены вершины графа, а потому в другой системе координат они попросту могут исчезнуть (но не появится ли одновременно других длинных дуг?). А с другой стороны, вне зависимости от системы координат их присутствие может быть сигналом о необходимости длительного хранения данных на определенном уровне иерархии, что накладывает дополнительные ограничения на эффективность реализации алгоритма. Одной из причин возникновения длинных дуг являются рассылки скалярных величин по всем итерациям какого-либо цикла: в таком виде длинные дуги не вызывают каких-либо серьезных проблем на практике.  
+
== Свойства алгоритма ==
 +
1. Один из самых быстрых алгоритмов вычисления вектора собственных значений и собственных векторов, соответствующих собственным значениям, матриц порядка большего, чем 25.
  
−
Для проектирования специализированных процессоров или реализации алгоритма на ПЛИС представляют интерес ''компактные укладки информационного графа'' [1], которые также имеет смысл привести в данном разделе.
+
2. Вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных – линейна.
 +
 +
3. Алгоритм не гарантирует высокой относительной точности для малых собственных значений.
  
 
= Программная реализация алгоритма =
 
= Программная реализация алгоритма =
−
Вторая часть описания алгоритмов в рамках AlgoWiki рассматривает все составные части процесса их реализации. Рассматривается как последовательная реализация алгоритма, так и параллельная. Описывается взаимосвязь свойств программ, реализующих алгоритм, и особенностей архитектуры компьютера, на которой они выполняются. Исследуется работа с памятью, локальность данных и вычислений, описывается масштабируемость и эффективность параллельных программ, производительность компьютеров, достигаемая на данной программе. Обсуждаются особенности реализации для разных классов архитектур компьютеров, приводятся ссылки на реализации в существующих библиотеках.
 
  
 
== Особенности реализации последовательного алгоритма ==
 
== Особенности реализации последовательного алгоритма ==
−
Здесь описываются особенности и варианты реализации алгоритма в виде последовательной программы, которые влияют на [[глоссарий#Эффективность реализации|''эффективность ее выполнения'']]. В частности, в данном разделе имеет смысл ''сказать о существовании блочных вариантов реализации алгоритма'', дополнительно описав потенциальные преимущества или недостатки, сопровождающие такую реализацию. Важный вопрос - это ''возможные варианты организации работы с данными'', варианты структур данных, наборов временных массивов и другие подобные вопросы. Для различных вариантов реализации следует оценить доступный ресурс параллелизма и объем требуемой памяти.
 
  
−
Важным нюансом является ''описание необходимой разрядности выполнения операций алгоритма'' (точности). На практике часто нет никакой необходимости выполнять все арифметические операции над вещественными числами с двойной точностью, т.к. это не влияет ни на устойчивость алгоритма, ни на точность получаемого результата. В таком случае, если значительную часть операций можно выполнять над типом float, и лишь в некоторых фрагментах необходим переход к типу double, это обязательно нужно отметить. Это прямое указание не только на правильную реализацию с точки зрения устойчивости по отношению к ошибкам округления, но и на более эффективную.
+
== Локальность данных и вычислений ==
  
−
Опираясь на информацию из [[#Описание ресурса параллелизма алгоритма|п.1.8]] (описание ресурса параллелизма алгоритма), при описании последовательной версии стоит сказать про возможности [[глоссарий#Эквивалентное преобразование|''эквивалентного преобразования программ'']], реализующих данных алгоритм. В дальнейшем, это даст возможность простого использования доступного параллелизма или же просто покажет, как использовать присущий алгоритму параллелизм на практике. Например, параллелизм на уровне итераций самого внутреннего цикла обычно используется для векторизации. Однако, в некоторых случаях этот параллелизм можно поднять "вверх" по структуре вложенности объемлющих циклов, что делает возможной и эффективную реализацию данного алгоритма на многоядерных SMP-компьютерах.
+
== Возможные способы и особенности параллельной реализации алгоритма ==
  
−
С этой же точки зрения, в данном разделе весьма полезны соображения по реализации алгоритма на различных параллельных вычислительных платформах. Высокопроизводительные кластеры, многоядерные узлы, возможности для векторизации или использования ускорителей - особенности этих архитектур не только опираются на разные свойства алгоритмов, но и по-разному должны быть выражены в программах, что также желательно описать в данном разделе.
+
== Масштабируемость алгоритма и его реализации ==
  
−
== Локальность данных и вычислений ==
+
Исследование проводилось на суперкомпьютере "Ломоносов" [http://parallel.ru/cluster Суперкомпьютерного комплекса Московского университета].
−
Вопросы локальности данных и вычислений не часто изучаются на практике, но именно локальность определяет эффективность выполнения программ на современных вычислительных платформах [2, 3]. В данном разделе приводятся оценки степени [[глоссарий#Локальность использования данных|''локальности данных'']] и [[глоссарий#Локальность вычислений|вычислений]] в программе, причем рассматривается как [[глоссарий#Временная локальность|''временна́я'']], так и [[глоссарий#Пространственная локальность|''пространственная'']] локальность. Отмечаются позитивные и негативные факты, связанные с локальностью, какие ситуации и при каких условиях могут возникать. Исследуется, как меняется локальность при переходе от последовательной реализации к параллельной. Выделяются ключевые шаблоны взаимодействия программы, реализующей описываемый алгоритм, с памятью. Отмечается возможная взаимосвязь между используемыми конструкциями языков программирования и степенью локальности, которыми обладают результирующие программы.  
+
 
 +
Для исследования масштабируемости алгоритм был реализован с использованием библиотеки OpenMP 1.5.5 для языка C++(компилятор g++) на cуперкомпьютере Ломоносов на узле compiler. Максимальное число потоков для данного узла равно 8. Для оценки масштабируемости, программа была запущена в 1, 2, 4, 8 потока. Код исследуемой реализации представлен по [https://github.com/Zero-ger/Tasks следующей ссылке]
  
−
Отдельно приводятся профили взаимодействия с памятью для вычислительных ядер и ключевых фрагментов. Если из-за большого числа обращений по общему профилю сложно понять реальную специфику взаимодействия программ с памятью, то проводится последовательная детализация и приводится серия профилей более мелкого масштаба.
+
Набор и границы значений изменяемых параметров запуска реализации алгоритма:
  
−
На рис.3 и рис.4 показаны профили обращения в память для программ, реализующих разложение Холецкого и быстрое преобразование Фурье, по которым хорошо видна разница свойств локальности у данных алгоритмов.
+
* число процессоров [1 : 8] (1,2,3,4,8);
 +
* размер матрицы [4 : 8192].
  
−
[[file:Cholesky_locality1.jpg|thumb|center|700px|Рис.3 Реализация метода Холецкого. Общий профиль обращений в память]]
+
Параметры компиляции и запуска:  
−
[[file:fft 1.PNG|thumb|center|700px|Рис.4 Нерекурсивная реализация БПФ для степеней двойки. Общий профиль обращений в память]]
 
  
−
== Возможные способы и особенности параллельной реализации алгоритма ==
+
ssh compiler
−
Раздел довольно обширный, в котором должны быть описаны основные факты и положения, формирующие параллельную программу. К их числу можно отнести:
+
module add slurm
−
* представленный иерархически ресурс параллелизма, опирающийся на структуру циклических конструкций и на граф вызовов программы;
+
module add openmpi/1.5.5-icc
−
* комбинацию (иерархию) массового параллелизма и параллелизма конечного;
+
g++ GetValues.cxx  -fopenmp -std=c++0x -o _scratch/mytask
−
* возможные способы распределения операций между процессами/нитями;
+
cd _scratch/
−
* возможные способы распределения данных;
+
sbatch -p test -n1 run ./mytask
−
* оценку количества операций, объёма и числа пересылок данных (как общего числа, так и в пересчёте на каждый параллельный процесс);
 
  
−
и другие.
+
Для управления количеством потоков использовалась операция OpenMP omp_set_num_threads();
  
−
В этом же разделе должны быть даны рекомендации или сделаны комментарии относительно реализации алгоритма с помощью различных технологий параллельного программирования: MPI, OpenMP, CUDA или использования директив векторизации.
 
  
−
== Масштабируемость алгоритма и его реализации ==
+
[[File:Graph.png|thumb|center|left|800px|Рис. 5. График времени работы программы в зависимости от размера задачи и числа процессов]]
−
Задача данного раздела - показать пределы [[глоссарий#Масштабируемость|''масштабируемости'']] алгоритма на различных платформах. Очень важный раздел. Нужно выделить, описать и оценить влияние точек барьерной синхронизации, глобальных операций, операций сборки/разборки данных, привести оценки или провести исследование [[глоссарий#Сильная масштабируемость|''сильной'']] и [[глоссарий#Слабая масштабируемость|''слабой'']] масштабируемости алгоритма и его реализаций.
 
  
−
Масштабируемость алгоритма определяет свойства самого алгоритма безотносительно конкретных особенностей используемого компьютера. Она показывает, насколько параллельные свойства алгоритма позволяют использовать возможности растущего числа процессорных элементов. Масштабируемость параллельных программ определяется как относительно конкретного компьютера, так и относительно используемой технологии программирования, и в этом случае она показывает, насколько может вырасти реальная производительность данного компьютера на данной программе, записанной с помощью данной технологии программирования, при использовании бóльших вычислительных ресурсов (ядер, процессоров, вычислительных узлов).
 
  
−
Ключевой момент данного раздела заключается в том, чтобы показать ''реальные параметры масштабируемости программы'' для данного алгоритма на различных вычислительных платформах в зависимости от числа процессоров и размера задачи  [4]. При этом важно подобрать такое соотношение между числом процессоров и размером задачи, чтобы отразить все характерные точки в поведении параллельной программы, в частности, достижение максимальной производительности, а также тонкие эффекты, возникающие, например, из-за блочной структуры алгоритма или иерархии памяти.
+
Проведем оценки масштабируемости:
  
−
На рис.5. показана масштабируемость классического алгоритма умножения плотных матриц в зависимости от числа процессоров и размера задачи. На графике хорошо видны области с большей производительностью, отвечающие уровням кэш-памяти.
+
* '''По размеру задачи''' — при увеличении размерности матрицы [[Глоссарий#Эффективность реализации|эффективность реализации]] постепенно уменьшается. Это связано с тем, что время чтения и подготовки исходных данных растёт быстрее, чем время работы исследуемого алгоритма.
−
[[file:Масштабируемость перемножения матриц производительность.png|thumb|center|700px|Рис.5 Масштабируемость классического алгоритма умножения плотных матриц в зависимости от числа процессоров и размера задачи]]
+
* '''По числу процессов''' - при увеличении числа процессов [[Глоссарий#Эффективность реализации|эффективность реализации]] увеличивается. Это связано с тем, что время чтения и подготовки исходной матрицы остаётся неизменным, а время работы алгоритма уменьшается, однако рост незначительный.
 +
* При одновременном увеличении числа процессов и размерности матрицы [[Глоссарий#Эффективность реализации|эффективность реализации]] уменьшается.
  
 
== Динамические характеристики и эффективность реализации алгоритма ==
 
== Динамические характеристики и эффективность реализации алгоритма ==
−
Это объемный раздел AlgoWiki, поскольку оценка эффективности реализации алгоритма требует комплексного подхода [5], предполагающего аккуратный анализ всех этапов от архитектуры компьютера до самого алгоритма. Основная задача данного раздела заключается в том, чтобы оценить степень эффективности параллельных программ, реализующих данный алгоритм на различных платформах, в зависимости от числа процессоров и размера задачи. Эффективность в данном разделе понимается широко: это и [[глоссарий#Эффективность распараллеливания|''эффективность распараллеливания'']] программы, это и [[глоссарий#Эффективность реализации|''эффективность реализации'']] программ по отношению к пиковым показателям работы вычислительных систем.
 
−
 
−
Помимо собственно показателей эффективности, нужно описать и все основные причины, из-за которых эффективность работы параллельной программы на конкретной вычислительной платформе не удается сделать выше. Это не самая простая задача, поскольку на данный момент нет общепринятой методики и соответствующего инструментария, с помощью которых подобный анализ можно было бы провести. Требуется оценить и описать эффективность работы с памятью (особенности профиля взаимодействия программы с памятью), эффективность использования заложенного в алгоритм ресурса параллелизма, эффективность использования коммуникационной сети (особенности коммуникационного профиля), эффективность операций ввода/вывода и т.п. Иногда достаточно интегральных характеристик по работе программы, в некоторых случаях полезно показать данные мониторинга нижнего уровня, например, по загрузке процессора, кэш-промахам, интенсивности использования сети Infiniband и т.п. Хорошее представление о работе параллельной MPI-программы дают данные трассировки, полученные, например, с помощью системы Scalasca.
 
  
 
== Выводы для классов архитектур ==
 
== Выводы для классов архитектур ==
−
В данный раздел должны быть включены рекомендации по реализации алгоритма для разных классов архитектур. Если архитектура какого-либо компьютера или платформы обладает специфическими особенностями, влияющими на эффективность реализации, то это здесь нужно отметить.
 
  
−
На практике это сделать можно по-разному: либо все свести в один текущий раздел, либо же соответствующие факты сразу включать в предшествующие разделы, где они обсуждаются и необходимы по смыслу. В некоторых случаях, имеет смысл делать отдельные варианты всей [[#ЧАСТЬ. Программная реализация алгоритмов|части II]] AlgoWiki применительно к отдельным классам архитектур, оставляя общей машинно-независимую [[#ЧАСТЬ. Свойства и структура алгоритмов|часть I]]. В любом случае, важно указать и позитивные, и негативные факты по отношению к конкретным классам. Можно говорить о возможных вариантах оптимизации или даже о "трюках" в написании программ, ориентированных на целевые классы архитектур.
+
== Существующие реализации алгоритма ==
  
−
== Существующие реализации алгоритма ==
+
Метод "Разделяй и властвуй" реализован в пакете LAPACK, описание которого можно найти в данной [http://www.netlib.org/lapack/lawnspdf/lawn132.pdf статье].
−
Для многих пар алгоритм+компьютер уже созданы хорошие реализации, которыми можно и нужно пользоваться на практике. Данный раздел предназначен для того, чтобы дать ссылки на основные существующие последовательные и параллельные реализации алгоритма, доступные для использования уже сейчас. Указывается, является ли реализация коммерческой или свободной, под какой лицензией распространяется, приводится местоположение дистрибутива и имеющихся описаний. Если есть информация об особенностях, достоинствах и/или недостатках различных реализаций, то это также нужно здесь указать. Хорошими примерами реализации многих алгоритмов являются MKL, ScaLAPACK, PETSc, FFTW, ATLAS, Magma и другие подобные библиотеки.
 
  
 
= Литература =
 
= Литература =
−
[1] Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. - СПб.: БХВ-Петербург, 2002. - 608 с.
+
[1] Джеймс Деммель. Вычислительная линейная алгебра. Теория и приложения. изд. Мир, 1/1/2001. — С. 228 - 236.
 +
 
 +
[2]  [https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_algorithm Divide-and-conquer eigenvalue algorithm]
  
−
[2] Воеводин В.В., Воеводин Вад.В. Спасительная локальность суперкомпьютеров //Открытые системы. - 2013. - № 9. - С. 12-15.
+
[3] Gr´egoire Pichon, Azzam Haidar, Mathieu Faverge, Jakub Kurzak. Divide and Conquer Symmetric Tridiagonal Eigensolver
 +
for Multicore Architectures. 2015. — С. 1 - 11.
  
−
[3] Воеводин Вад.В., Швец П. Метод покрытий для оценки локальности использования данных в программах // Вестник УГАТУ. — 2014. — Т. 18, № 1(62). — С. 224–229.
+
[4] C.F. Borges, W.B. Gragg, A parallel divide-and-conquer method for the generalized real symmetric definite tridiagonal
 +
eigenproblem, in: L. Reichel, A. Ruttan, R.S. Varga (Eds.), Numerical Linear Algebra, Proc. Conf. in Numerical
 +
Linear Algebra and Scientific Computation, Kent, OH, de Gruyter, Berlin, 1993. — С. 11–29
  
−
[4] Антонов А.С., Теплов А.М. О практической сложности понятия масштабируемости параллельных программ// Высокопроизводительные параллельные вычисления на кластерных системах (HPC 2014): Материалы XIV Международной конференции -Пермь: Издательство ПНИПУ, 2014. С. 20-27.
+
[5] L. Elsner, A. Fasse, E. Langmann. A divide-and-conquer method for the tridiagonal generalized eigenvalue problem. Journal of Computational and Applied Mathematics. 1997. — С. 141 - 148.
  
−
[5] Никитенко Д.А. Комплексный анализ производительности суперкомпьютерных систем, основанный на данных системного мониторинга // Вычислительные методы и программирование. 2014. 15. 85–97.
+
[6] [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]
  
 
[[en:Description of algorithm properties and structure]]
 
[[en:Description of algorithm properties and structure]]

Текущая версия на 17:50, 3 февраля 2017

Symbol wait.svgЭта работа прошла предварительную проверку
Дата последней правки страницы: 03.02.2017
Данная работа соответствует формальным критериям.
Проверено ASA.


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

Авторы:

Изотова И.А. (619 группа) пункты 1.1, 1.2, 1.5, 1.7, 1.9

Ульянов Г.В. (619 группа) пункты 1.3, 1.4, 1.6, 1.8, 1.10


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

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

«Разделяй и властвуй» - метод вычисления собственных значений и векторов симметричной трехдиагональной матрицы. В настоящее время это один из самых быстрых методов при вычислении собственных значений и векторов симметричной трехдиагональной матрицы порядка [math]n \gt 25[/math] (точное значение этого порогового порядка зависит от компьютера.) Однако алгоритм не гарантирует высокой относительной точности для малых собственных значений.Его численно устойчивая реализация весьма не тривиальна. Хотя впервые метод был предложен еще в 1981 г., "правильный" способ его реализации был найден лишь в 1992 г. Этот способ воплощен LAPACK-программами ssyevd (для плотных матриц) и sstevd (для трехдиагональных матриц). В них стратегия "разделяй-и-влавствуй" используется для матриц порядка, большего чем 25. Для матриц меньшего порядка (или если нужны только собственные значения) происходит автоматический переход к QR-итерации. Однако, из-за высокой сложности, алгоритм пока не вошел в состав ALGLIB.

В наихудшем случае алгоритм «разделяй-и-властвуй» требует [math] O(n^{3})[/math] операций, но константа, на практике оказывается весьма малой. Для некоторых специальных распределений собственных значений — [math] O(n^{2})[/math] (подробнее в разделе 1.6 Последовательная сложность алгоритма).

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

Разделяй

Пусть матрицы [math]T, T_{1}[/math], [math]T_{2}[/math] имеют размерность [math]n \times n[/math], [math]m \times m[/math], [math](n - m) \times (n - m)[/math].

Запишем матрицу [math]T[/math] в следующем виде:

[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]\hat{T}_{1}[/math] и [math]\hat{T}_{2}[/math], т.е.: [math]\hat{T}_{1} = Q_{1} D_{1} Q_{1}^{T}[/math] и [math]\hat{T}_{2} = Q_{2} D_{2} Q_{2}^{T}[/math].

Властвуй

Установим связь между собственными значениями матрицы [math]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}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] Q_{1}^{T}[/math], а во второй - элементы первого столбца матрицы [math] Q_{2}^{T}[/math], так как [math]v = \begin{bmatrix} 0 , \ldots , 0 , 1 , 1 , 0 \ldots , 0 \end{bmatrix}^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} \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]

По Рис.1 можно заметить, что прямая [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 = 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 = \alpha [(D - \alpha I)^{-1}u][/math],

поскольку

[math] \rho u^{T}(D - \alpha I)^{-1}u + 1 = f(\alpha) = 0 [/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\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] достаточно мало.

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] о чём будет говориться ниже, в разделе 1.6 Последовательная сложность алгоритма.

Данной операции предшествует вычисление собственных значений и векторов матрицы [math] D + \rho uu^{T}[/math].

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

Алгоритм "разделяй и властвуй" для вычисления собственных значений и векторов симметричной трехдиагональной матрицы состоит из вызова рекурсивной функции proc dc_eig [math](T, Q, \Lambda)[/math], выполнение которой описывается следующими этапами:

  • Подаваемая на вход матрица [math]T[/math] представляется в виде [math] T = \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T} [/math], представление в таком виде возможно до тех пор, пока размерность матрицы [math]T[/math] не достгнет 1.
  • От матриц [math]T_{1}[/math] и [math]T_{2}[/math] также вызываются соотвтетствующие функции proc dc_eig [math](T_{1}, Q_{1}, \Lambda_{1})[/math] и proc dc_eig [math](T_{2}, Q_{2}, \Lambda_{2})[/math].
  • Вычисляются матрицы [math]D+\rho uu^{T}[/math] по [math] \Lambda_{1},\Lambda_{2}, Q_{1}, Q_{2}[/math].
  • Вычисление собственных значений[math]\Lambda[/math] и векторов [math]Q^{'}[/math] для матрицы [math]D+\rho uu^{T}[/math].
  • Построение матрицы собственных векторов [math]Q[/math]для матрицы [math]T[/math].

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

Рассмотри реализацию метода "разделяй и властвуй" в виде алгоритма, состоящего из частного случая, когда матрица [math]T[/math] имеет размерность [math]1 \times 1[/math] и общего случая, когда матрица [math]T[/math] - симметричная трехдиагональная матрица произвольной размерности, кроме [math]1 \times 1[/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]
    if матрица [math]T[/math] имеет размерность [math]1 \times 1[/math]                       (проверяем матрицу [math]T[/math] на принадлежность частному случаю алгоритма)
    then
        присвоить выходным параметрам значения [math]Q = 1, \Lambda = T[/math]              (если матрица [math]T[/math] [math]1 \times 1[/math], то СЗ и СВ матрицы [math]T[/math] найдены)
    else                                                                                    (если матрица [math]T[/math] имеет общей вид, то:)
        [math]T[/math] в виде [math] T = \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T} [/math]                                                                  (разложим матрицу [math]T[/math])
        call dc_eig[math](T_{1},Q_{1},\Lambda_{1})[/math]                              (вызываем рекурсивно функцию dc_eig, вычисляющую СЗ и СВ матрицы [math]T_1[/math])
        call dc_eig[math](T_{2},Q_{2},\Lambda_{2})[/math]                              (вызываем рекурсивно функцию dc_eig, вычисляющую СЗ и СВ матрицы [math]T_2[/math])
        Построить [math]D+\rho uu^{T}[/math] по [math] \Lambda_{1},\Lambda_{2}, Q_{1}, Q_{2}[/math]
        Найти матрицу собственных значений [math]\Lambda[/math]
        Найти матрицу собственных векторов [math]Q^{'}[/math] для матрицы [math]D+\rho uu^{T}[/math]
        Построить матрицу собственных векторов [math]Q[/math] для матрицы [math]T[/math] : [math] Q = \begin{bmatrix} Q_{1} & 0 \\ 0 & 
Q_{2}\end{bmatrix}* Q^{'} [/math]
        Присвоить выходным параметрам значения [math]Q[/math] и [math]\Lambda[/math]
 endif

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

Пусть [math]t(n)[/math] - число флопов при обработке матрицы размера [math]n \times 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\lt \lt 1[/math] ( так как матрица [math]Q^{'}[/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 = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}*Q^{'}[/math] является самым затратным действием, то дефляция многих собственных значений может привести в снижению сложности до [math]O(n^{2})[/math] операций.

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

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

Рис. 2 Информационный граф шага "разделяй"
Рис. 3 Связь шагов "Разделяй" и "властвуй"
Рис. 4 Информационный граф шага "властвуй"

Операции:

  • [math]D[/math] - операция "разделяй" разделения матрицы [math] T [/math] на матрицы [math]T_{1}[/math] и [math]T_{2}[/math]
  • [math]C_{n}[/math] - операция "властвуй" сбора матрицы [math] D_{n*n} [/math]
  • [math] R_{n} [/math] - операция поиска собственных значений матрицы [math] T_{n*n} [/math]
  • [math] S^{'}_{n} [/math] - операция поиска собственных векторов матрицы [math] D_{n*n} [/math]
  • [math] S_{n} [/math] - операция поиска собственных векторов матрицы [math] T_{n*n} [/math]

Данные:

  • [math]Q_{n}[/math] - вектор собственных значений размерности n.
  • [math]V_{n*n}[/math] - матрица собственных векторов размерности n.

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

Параллелизм достигается за счёт вызова рекурсивных функций dc_eig, вычисляющих собственные значения и собственные вектора матриц [math]T_1[/math] и [math]T_2[/math]. Рекурсивное разделение матрицы приводит к возможности построения структуры в виде двоичного дерева, в каждой вершине которого происходит вычисение функции dc_eig. Таким образом высота ярусно-параллельной формы, где вычисления для каждого листа дерева в такой структуре происходят независимо, составляет [math]log_2(n)[/math].

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

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

Матрица [math]T[/math] размерностью [math]n \times n[/math]:

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

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

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

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

Собственные значения [math]\alpha_{i}[/math] (всего [math]n[/math]) и собственные векторы [math]\lambda_{i}[/math] (всего [math]n^2[/math]).

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

[math]n+n^2[/math]

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

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

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

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

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

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

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

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

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

Исследование проводилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского университета.

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

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

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

Параметры компиляции и запуска: 

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

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


Рис. 5. График времени работы программы в зависимости от размера задачи и числа процессов


Проведем оценки масштабируемости:

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

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

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

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

Метод "Разделяй и властвуй" реализован в пакете LAPACK, описание которого можно найти в данной статье.

3 Литература

[1] Джеймс Деммель. Вычислительная линейная алгебра. Теория и приложения. изд. Мир, 1/1/2001. — С. 228 - 236.

[2] Divide-and-conquer eigenvalue algorithm

[3] Gr´egoire Pichon, Azzam Haidar, Mathieu Faverge, Jakub Kurzak. Divide and Conquer Symmetric Tridiagonal Eigensolver for Multicore Architectures. 2015. — С. 1 - 11.

[4] C.F. Borges, W.B. Gragg, A parallel divide-and-conquer method for the generalized real symmetric definite tridiagonal eigenproblem, in: L. Reichel, A. Ruttan, R.S. Varga (Eds.), Numerical Linear Algebra, Proc. Conf. in Numerical Linear Algebra and Scientific Computation, Kent, OH, de Gruyter, Berlin, 1993. — С. 11–29

[5] L. Elsner, A. Fasse, E. Langmann. A divide-and-conquer method for the tridiagonal generalized eigenvalue problem. Journal of Computational and Applied Mathematics. 1997. — С. 141 - 148.

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

Источник — https://algowiki-project.org/w/ru/index.php?title=Участник:Amirida/Метод_«разделяй_и_властвуй»_вычисления_собственных_значений_и_векторов_симметричной_трехдиагональной_матрицы&oldid=22103