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

Файл:Symbol wait.svg

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

Авторы статьи: Сергеев Никита (группа 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)

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

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

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

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

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

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

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

Пусть

[math]\displaystyle{ L = \begin{bmatrix} a_{1} & b_{1}&&&&& \\ b_{1} & \ddots & \ddots \\ & \ddots & a_{m-1} & b_{m-1} \\ && b_{m-1} & a_{m} & b_{m} \\ &&& b_{m} & a_{m+1} & b_{m+1} \\ &&&& b_{m+1} & \ddots \\ &&&&&& \ddots & b_{n-1} \\ &&&&&& b_{n-1} & a_{n} \\ \end{bmatrix} = \begin{bmatrix} a_{1} & b_{1} &&&&& \\ b_{1} & \ddots & \ddots \\ & \ddots & a_{m-1} & b_{m-1} \\ && b_{m-1} & a_{m} - b_{m} \\ &&&& a_{m+1} - b_{m} & b_{m+1} \\ &&&& b_{m+1} & \ddots \\ &&&&&& \ddots & b_{n-1} \\ &&&&&& b_{n-1} & a_{n} \\ \end{bmatrix} + }[/math]

[math]\displaystyle{ + \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]\displaystyle{ T_{1} }[/math] и [math]\displaystyle{ T_{2} }[/math]: [math]\displaystyle{ T_{i}=Q_{i} \Lambda_{i} Q_{i}^{T} }[/math]. На самом деле, они будут вычисляться рекурсивно тем же самым алгоритмом.

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

[math]\displaystyle{ T = \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T} = \begin{bmatrix} Q_{1} \Lambda_{1} Q_{1}^{T} & 0 \\ 0 & Q_{2} L_{2} Q_{2}^{T}\end{bmatrix} + b_{m}vv^{T} = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}(\begin{bmatrix} \Lambda_{1} & \\ & \Lambda_{2}\end{bmatrix} + b_{m}vv^{T})\begin{bmatrix} Q_{1}^{T} & 0 \\ 0 & Q_{2}^{T}\end{bmatrix} }[/math],

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

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

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

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

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

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

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

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

 Лемма 2. Если [math]\displaystyle{ \alpha }[/math] - собственное значение матрицы [math]\displaystyle{ D + \rho uu^{T} }[/math], то соответствующий вектор равен [math]\displaystyle{ (D - \alpha I)^{-1}u }[/math]. Поскольку матрица [math]\displaystyle{ D - \alpha I }[/math] диагональная, для вычисления такого вектора достаточно [math]\displaystyle{ O(n) }[/math] флопов.

Доказательство. [math]\displaystyle{ (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]\displaystyle{ = u + \alpha(D - \alpha I)^{-1}u - u = \alpha [(D - \alpha I)^{-1}u] }[/math],

поскольку

[math]\displaystyle{ \rho u^{T}(D - \alpha I)^{-1}u + 1 = f(\alpha) = 0 }[/math].

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

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

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

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

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

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

Предположим, что некоторое [math]\displaystyle{ u_{i} }[/math], хотя и мало, все же недостаточно мало для того, чтобы была зарегистрирована дефляция. В этом случае применение метода Ньютона к решению векового уравнения встречается с затруднениями. Вспомним, что пересчет приближённого решения [math]\displaystyle{ u_{j} }[/math] уравнения [math]\displaystyle{ f(\lambda) = 0 }[/math] в методе Ньютона основан на следующих положениях:

1. Вблизи точки [math]\displaystyle{ \lambda = \lambda_{j} }[/math] функция [math]\displaystyle{ f(\lambda) }[/math] аппроксимируется линейной функцией [math]\displaystyle{ l(\lambda) }[/math]; график есть прямая линия, касающаяся графика функции [math]\displaystyle{ f(\lambda) }[/math] при [math]\displaystyle{ \lambda = \lambda_{j} }[/math].

2. В качестве [math]\displaystyle{ \lambda_{j+1} }[/math] берётся нуль этого линейного приближения, т.е. [math]\displaystyle{ l(\lambda_{j+1})=0 }[/math].

Функция, показанная на рис.1, не доставляет видимых трудностей методу Ньютона, поскольку вблизи каждого своего нуля [math]\displaystyle{ f(\lambda) }[/math] достаточно хорошо аппроксимируется линейными функциями. Однако рассмотрим график функции на рис. 2. Она получена из функции на рис. 1 заменой значения .5 для [math]\displaystyle{ u_{i}^{2} }[/math] на .001. Это новое значение недостаточно мало для того, чтобы вызвать дефляцию. График функции в левой части рис.2 визуально не отличим от её вертикальных и горизонтальных асимптот, поэтому в правой части укрупненно воспроизведён фрагмент графика, прилегающий к вертикальной асимптоте [math]\displaystyle{ \lambda = 2 }[/math]. Видно, что график слишком быстро "выполняет поворот" и для большей части значений [math]\displaystyle{ \lambda }[/math] почти горизонтален. Поэтому, применяя метод Ньютона почти к любому начальному приближению [math]\displaystyle{ \lambda_{0} }[/math], мы получаем линейное приближение [math]\displaystyle{ l(\lambda) }[/math] с почти горизонтальным графиком и малым положительным угловым коэффициентом. В результате [math]\displaystyle{ \lambda_{1} }[/math] является отрицательным числом, огромным по абсолютной величине, которое совершенно бесполезно в качестве приближения к истинному корню.

Рис. 2. График функции [math]\displaystyle{ 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]\displaystyle{ f(\lambda) }[/math] нельзя хорошо приблизить линейной функцией [math]\displaystyle{ l(\lambda) }[/math], попробуем взять в качестве приближения какую-нибудь другую простую функцию [math]\displaystyle{ h(\lambda) }[/math]. Нет ничего особого именно в прямых линиях: для метода Ньютона вместо [math]\displaystyle{ l(\lambda) }[/math] можно взять любое приближение [math]\displaystyle{ h(\lambda) }[/math], значения и нули которого легко вычисляются. Функция [math]\displaystyle{ f(\lambda) }[/math] имеет полюсы в точках [math]\displaystyle{ d_{i} }[/math] и [math]\displaystyle{ d_{i+1} }[/math], которые определяют её поведение в соответствующих окрестностях. Поэтому при поиске корня в интервале [math]\displaystyle{ (d_{i+1}, d_{i}) }[/math] естественно выбрать функцию [math]\displaystyle{ h(\lambda) }[/math], также имеющую эти полюсы, т.е. функцию вида [math]\displaystyle{ h(\lambda)= \frac{c_{1}}{d_{i}-\lambda} + \frac{c_{2}}{d_{i+1}-\lambda} + c_{3} }[/math]

Константы [math]\displaystyle{ c_{1},c_{2} }[/math] и [math]\displaystyle{ c_{3} }[/math] обеспечивающие, что [math]\displaystyle{ h(\lambda) }[/math] есть приближение к [math]\displaystyle{ f(\lambda) }[/math], можно выбрать несколькими способами. Отметим, что если [math]\displaystyle{ c_{1},c_{2} }[/math] и [math]\displaystyle{ c_{3} }[/math] уже известны, то уравнение [math]\displaystyle{ h(\lambda)=0 }[/math] легко решается относительно [math]\displaystyle{ \lambda }[/math], поскольку сводится к эквивалентному квадратному уравнению [math]\displaystyle{ c_{1}(d_{i+1}-\lambda)+c_{2}(d_{i}-\lambda)+c_{3}(d_{i}-\lambda)(d_{i+1}-\lambda)=0 }[/math]

Пусть [math]\displaystyle{ \lambda_{j} }[/math] - приближённое значение корня. определим [math]\displaystyle{ c_{1},c_{2} }[/math] и [math]\displaystyle{ c_{3} }[/math] так, чтобы [math]\displaystyle{ \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]\displaystyle{ \lambda }[/math] в окрестности [math]\displaystyle{ \lambda_{j} }[/math]. Заметим, что [math]\displaystyle{ 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]\displaystyle{ \lambda \in (d_{i+1},d_{i}) }[/math], то [math]\displaystyle{ \psi_{1}(\lambda) }[/math] есть сумма положительных слагаемых, а [math]\displaystyle{ \psi_{2}(\lambda) }[/math] - сумма отрицательных. Поэтому и [math]\displaystyle{ \psi_{1}(\lambda) }[/math], и [math]\displaystyle{ \psi_{2}(\lambda) }[/math] могут быть вычислены с высокой точностью; однако при их сложении вполне вероятно взаимное уничтожение верных разрядов и потеря относительной точности в сумме. Возьмем числа [math]\displaystyle{ c_{1} }[/math] и [math]\displaystyle{ \hat{c_{1}} }[/math], такие, что функция [math]\displaystyle{ h_{1}(\lambda) \equiv \hat{c_{1}} + \frac{c_{1}}{d_{i}-\lambda} }[/math] удовлетворяет условиям [math]\displaystyle{ h_{1}(\lambda_{j}) = \psi_{1}(\lambda_{j}) }[/math] и [math]\displaystyle{ h_{1}^{'}(\lambda_{j})=\psi_{1}^{'}(\lambda_{j}) }[/math] (*)

Это означает, что гипербола, являющаяся графиком функции [math]\displaystyle{ h_{1}(\lambda) }[/math], касается графика функции [math]\displaystyle{ \psi_{i}(\lambda) }[/math] при [math]\displaystyle{ \lambda = \lambda_{j} }[/math]. Два условия в (*) - это обычные условия метода Ньютона, за исключением того, что вместо прямой в качестве приближения используется гипербола. Легко проверить, что [math]\displaystyle{ c_{1}=\psi_{1}^{'}(\lambda_{j})(d_{i} - \lambda_{j})^{2} }[/math] и [math]\displaystyle{ \hat{c_{1}}=\psi_{1}(\lambda_{j}) - \psi_{1}^{'}(\lambda_{j})(d_{i} - \lambda_{j}) }[/math].

Подобным же образом выбираем [math]\displaystyle{ c_{2} }[/math] и [math]\displaystyle{ \hat{c_{2}} }[/math] так, чтобы функция [math]\displaystyle{ h_{2}(\lambda) \equiv \hat{c_{2}} + \frac{c_{2}}{d_{i+1}-\lambda} }[/math] удовлетворяла условиям [math]\displaystyle{ h_{2}(\lambda_{j}) = \psi_{2}(\lambda_{j}) }[/math] и [math]\displaystyle{ h_{2}^{'}(\lambda_{j})=\psi_{2}^{'}(\lambda_{j}) }[/math]

Наконец, полагаем [math]\displaystyle{ 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].

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

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

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

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

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

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

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

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

В случае, если [math]\displaystyle{ T }[/math] - матрица размера [math]\displaystyle{ 1 }[/math] x [math]\displaystyle{ 1 }[/math] :

1. Присвоить выходным параметрам значения [math]\displaystyle{ Q = 1, \Lambda = T }[/math]

Иначе:

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

2. call dc_eig[math]\displaystyle{ (T_{1},Q_{1},\Lambda_{1}) }[/math]

3. call dc_eig[math]\displaystyle{ (T_{2},Q_{2},\Lambda_{2}) }[/math]

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

5. Найти матрицу собственных значений [math]\displaystyle{ \Lambda }[/math]

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

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

8. Присвоить выходным параметрам значения [math]\displaystyle{ Q }[/math] и [math]\displaystyle{ \Lambda }[/math]

endif

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

Рассмотрим сложность алгоритма. Пусть [math]\displaystyle{ t(n) }[/math] - число флопов, которое потребуется при обработке матрицы размера [math]\displaystyle{ n * n }[/math] процедурой dc_eig. В таком случае

 [math]\displaystyle{  t(n) = 2t(n/2)  }[/math] два рекурсивных обращения к dc_eig[math]\displaystyle{ (T_{i},Q_{i},\Lambda_{i}) }[/math]
        [math]\displaystyle{  +O(n^{2}) }[/math] вычисление собственных значений матрицы [math]\displaystyle{ D+\rho uu^{T} }[/math]
        [math]\displaystyle{  +O(n^{2}) }[/math] вычисление собственных векторов матрицы [math]\displaystyle{ D+\rho uu^{T} }[/math]
        [math]\displaystyle{  +c*n^{3} }[/math] вычисление матрицы [math]\displaystyle{ Q = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}*Q^{'} }[/math]

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

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

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

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

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

Операции:

[math]\displaystyle{ D }[/math] - операция Divide разделения матрицы [math]\displaystyle{ T }[/math] на матрицы [math]\displaystyle{ T_{1} }[/math] и [math]\displaystyle{ T_{2} }[/math]

Индексы - размерность

[math]\displaystyle{ S_{k} }[/math] - операция Conquer сбора матрицы [math]\displaystyle{ D }[/math]

[math]\displaystyle{ L_{k} }[/math] - операция поиска собственных значений

[math]\displaystyle{ S^{'} }[/math] - операция поиска собственных векторов матрицы [math]\displaystyle{ Q^{'} }[/math]

[math]\displaystyle{ S }[/math] - операция поиска собственных векторов матрицы [math]\displaystyle{ Q }[/math]

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

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

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

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

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

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

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

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

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

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

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

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

[math]\displaystyle{ L = \begin{bmatrix} a_{1} & b_{1}&&&&& \\ b_{1} & \ddots & \ddots \\ & \ddots & a_{m-1} & b_{m-1} \\ && b_{m-1} & a_{m} & b_{m} \\ &&& b_{m} & a_{m+1} & b_{m+1} \\ &&&& b_{m+1} & \ddots \\ &&&&&& \ddots & b_{n-1} \\ &&&&&& b_{n-1} & a_{n} \\ \end{bmatrix} }[/math]

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

Выходные данные: Вектор собственных значений длины [math]\displaystyle{ n }[/math] и [math]\displaystyle{ n }[/math] собственных векторов длины [math]\displaystyle{ n }[/math] , соответствующих собственным значениям исходной матрицы.

Объём выходных данных: [math]\displaystyle{ n(n + 1) }[/math]

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

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

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

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

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

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

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

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

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

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

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

Тестирование реализации алгоритма проводилось на суперкомпьютере "Ломоносов" с использованием 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

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

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

Литература

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

[2] Doron Grill and Eitan Tadmor AN [math]\displaystyle{ O(N2) }[/math] METHOD FOR COMPUTING THE EIGENSYSTEM OF [math]\displaystyle{ N }[/math]x[math]\displaystyle{ N }[/math] 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