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

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

Авторы:

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

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 Информационный граф шага "властвуй"

Операции:

• [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.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();

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