Участник:Vika3410/Метод Разделяй и властвуй

Эта работа ждет рассмотрения преподавателем Дата последней правки страницы: 31.12.2016 Авторы этой статьи считают, что задание выполнено.

Авторы статьи: Сергеев Никита (группа 614) (1.1, 1.2, 1.5, 1.6, 1.7, 1.8, 1.10, 2.4), Стрельцова Виктория (группа 614) (1.2, 1.3, 1.4, 1.6, 1.7, 1.9, 1.10, 2.4)

Метод "Разделяй и властвуй" вычисления собственных значений и векторов трёхдиагональной матрицы
Последовательный алгоритм
Последовательная сложность	$c\frac{4}{3}n^{3}$
Объём входных данных	$2n + 1$
Объём выходных данных	$n(n + 1)$
Параллельный алгоритм
Параллельная сложность	$O(n^{2.3})$ крайне редко: $O(n^{2})$

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

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

Метод разделяй и властвуй вычисления собственных значений и векторов трёхдиагональной матрицы - является наиболее быстрым из существующих методов вычисления всех собственных значений и собственных векторов трехдиагональной матрицы, начиная с порядка n, который примерно равен 26. (Точное значение этого порогового порядка в конкретном случае зависит от компьютера.)

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

Для матриц порядка ниже 25 (или в тех случаях, когда требуются только собственные значения) автоматически происходит переход к методу QR-итерации.

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

Пусть

$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} +$

$+ \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}$

В предположении, что известны спектральные разложения матриц $T_{1}$ и $T_{2}$: $T_{i}=Q_{i} \Lambda_{i} Q_{i}^{T}$. На самом деле, они будут вычисляться рекурсивно тем же самым алгоритмом.

Можно установить связь между собственными значениями матрицы Т и собственными значениями матриц $T_{1}$ и $T_{2}$. Получается:

$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}$,

где $u = \begin{bmatrix} Q_{1}^{T} & 0 \\ 0 & Q_{2}^{T}\end{bmatrix}v$ - матрица, в первой строке которой находятся элементы последнего столбца матрицы $Q_{1}^{T}$, а во второй - элементы первого столбца матрицы $Q_{2}^{T}$, так как $v = \begin{bmatrix} 0 , \ldots , 0 , 1 , 1 , 0 \ldots , 0 \end{bmatrix}^T$. Следовательно, $T$ имеет такие же собственные значения, что и подобная ей матрица $D + \rho uu^{T}$, где $D = \begin{bmatrix} L_{1} & 0 \\ 0 & L_{2}\end{bmatrix}$ - диагональная матрица, $\rho = b_{m}$ - число, а $u$ - вектор. В дальнейшем предполагаем, не ограничивая общности, что диагональные элементы $d_{1}, \ldots, d_{n}$ матрицы $D$ упорядочены по не возрастанию: $d_{n} \lt = \ldots \lt =d_{1}$.

Для того, чтобы найти собственные значения матрицы $D + \rho uu^{T}$, потребуется вычислить её характеристический многочлен, считая пока матрицу $D - \lambda I$ невырожденной. Тогда получаем $det(D + \rho uu^{T} - \lambda I) = det((D - \lambda I)(I + \rho (D- \lambda I)^{-1} uu^{T}))$.

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

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

Следовательно, получаем $det(I + \rho (D - \lambda I)^{-1}uu^{T}) = 1 + \rho u^{T}(D - \lambda I)^{-1}u$ $= 1 + \rho \sum_{i=1, n} \frac{u_{i}^{2}} {d_{i}-\lambda} \equiv f(\lambda)$ ,

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

Рис. 1. График функции $f(\lambda) = 1 + \frac{0.5}{1 - \lambda} + \frac{0.5}{2 - \lambda} + \frac{0.5}{3 - \lambda} + \frac{0.5}{4 - \lambda}$

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

 Лемма 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] флопов.

Доказательство. $(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]=$

$= u + \alpha(D - \alpha I)^{-1}u - u = \alpha [(D - \alpha I)^{-1}u]$,

поскольку

$\rho u^{T}(D - \alpha I)^{-1}u + 1 = f(\alpha) = 0$.

Что и требовалось доказать.

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

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

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

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

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

Вычислительным ядром последовательной схемы решения будет являться вычисление матрицы $Q$ собственных векторов при помощи умножения матрицы $Q = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}$ на матрицу $Q^{'}$. Эта операция имеет сложность $cn^{3}$ о чём будет говориться ниже, в разделе #Последовательная сложность алгоритма .

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

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

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

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

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

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

В случае, если $T$ - матрица размера $1$ x $1$ :

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

Иначе:

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

2. call dc_eig$(T_{1},Q_{1},\Lambda_{1})$

3. call dc_eig$(T_{2},Q_{2},\Lambda_{2})$

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

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

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

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

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

endif

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

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

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

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

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

На рисунке 3 изображена структура алгоритма.

Рис. 3. Дерево алгоритма "Разделяй и властвуй"

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

В структуре алгоритма имеем два параллельных блока - вызов рекурсивных функций dc_eig для вычисления собственных значений и векторов матриц $T_{1}$ и $T_{2}$ - это единственная часть алгоритма, в которой возможно использовать параллелизм.

Рекурсивное разделение матрицы ведёт к иерархии вычислительных подблоков с графом зависимости данных, который принимает вид двоичного дерева. Такая структура является естественным путем разделения вычислений между процессорами.

- В вершине дерева все процессы взаимодействуют.

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

- В листьях деревьев каждый процесс занимается вычислением независимо.

В случае реализации без параллелизма, функции dc_eig отрабатывают последовательно - сначала со входными параметрами $T_{1}, Q_{1}, \Lambda_{1}$, затем - $T_{2}, Q_{2}, \Lambda_{2}$.

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

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

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

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

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

$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}$

Объём входных данных: $2n + 1$ (в силу симметричности достаточно хранить только диагональ и над/под диагональные элементы).

Выходные данные: Вектор собственных значений длины $n$ и $n$ собственных векторов длины $n$ , соответствующих собственным значениям исходной матрицы.

Объём выходных данных: $n(n + 1)$

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

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

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

Она удовлетворяет соотношению $t(n) = n^3 +2t(\frac{n}{2})$

Параллельная сложность алгоритма (за счет дефляции): $O(n^{2.3})$; в крайне редких случаях: $O(n^{2})$.

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

1. Дефляция (#Дефляция)

2. Адаптированный метод Ньютона (#Макроструктура алгоритма)

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

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

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

Тестирование реализации алгоритма проводилось на суперкомпьютере "Ломоносов" с использованием 1, 2, 4, 8, 16 процессоров на матрицах размера от 4 до 4096. Реализация произведена на библиотеке OpenMP.

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

Ниже приведена схема зависимости времени работы программы от размера матрицы при разном числе процессоров. При размере матрицы от 4 до 1024 время работы очень мало, поэтому эти данные не указаны на схеме.

Рис. 4. Зависимость времени работы от размера матрицы

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

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

Код реализации: https://drive.google.com/open?id=0B0G2UV7cuLwHQk9xd1lyeDNJOXc

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

Единственный найденный пример реализации описан в Francoise Tisseury, Jack Dongarra, "Parallelizing the Divide and Conquer Algorithm for the Symmetric Tridiagonal Eigenvalue Problem on Distributed Memory Architectures" ([3] с. 7, 19-20).

3 Литература

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

[2] Doron Grill and Eitan Tadmor AN $O(N2)$ METHOD FOR COMPUTING THE EIGENSYSTEM OF $N$x$N$ SYMMETRIC TRIDIAGONAL MATRICES BY THE DIVIDE AND CONQUER APPROACH

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

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

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