Участник:Vika3410/Метод Разделяй и властвуй
1 ЧАСТЬ. Свойства и структура алгоритмов[править] 1.1 Общее описание алгоритма[править] Метод разделяй и властвуй вычисления собственных значений и векторов трёхдиагональной матрицы - это наиболее быстрый из существующих методов, если нужны все собственные значения и собственные векторы трехдиагональной матриц, начиная с порядка n, примерно равного 26. (Точное значение этого порогового порядка зависит от компьютера.) Его численно устойчивая реализация весьма не тривиальна. В самом деле, хотя впервые метод был предложен еще в 1981 г., "правильный" способ его реализации был найден лишь в 1992 г. Этот способ воплощен LAPACK-программами ssyevd (для плотных матриц) и sstevd (для трехдиагональных матриц). В них стратегия "разделяй-и-влавствуй" используется для матриц порядка, большего чем 25. Для матриц меньшего порядка (или если нужны только собственные значения) происходит автоматический переход к QR-итерации. 1.2 Математическое описание алгоритма[править] Пусть
Предположим, что нам известны спектральные разложения матриц и : . В действительности, они будут рекурсивно вычисляться тем же самым алгоритмом. Установим связь между собственными значениями матрицы Т и собственными значениями матриц и . Имеем:
,
где
так как , получим матрицу, состоящую из последнего столбца матрицы и первого столбца матрицы . Следовательно, имеет те же собственные значения, что и пдобная ей матрица , где - диагональная матрица, - число, а - вектор. Будем предполагать, не ограничивая общности, что диагональные элементы матрицы упорядочены по убыванию: . Чтобы найти собственные значения матрицы , вычислим её характеристический многочлен, считая пока матрицу невырожденной. Тогда
.
Поскольку невырожденна, тогда и только тогда, когда - собственное значение. Заметим, что матрица получается из единичной добавлением матрицы ранга 1. Определитель такой матрицы легко вычислить.
Лемма 1. Справедливо равенство , где и - векторы.
Таким образом,
,
т.е. собственные значения матрицы есть корни так называемого векового уравнения . Если все числа различны и все (случай общего положения), то имеет график типа показанного на рис.1(где и ).
Рис. 1. График функции Можно видеть, что прямая является горизонтальной асимптотой для этого графика, а прямые есть вертикальные асимптоты. Поскольку , функция возрастает всюду, кроме точек . Поэтому корни функции разделяются числами и ещё один корень находится справа от точки (на рис. 1 ). (При функция всюду убывает и соответствующий корень находится слева от точки ). Для функции , монотонной и гладкой на каждом из интервалов , можно найти вариант метода Ньютона, который быстро и монотонно сходится к каждому из корней, если начальная точка взята в . Нам достаточно знать, что на практике метод сходится к каждому собственному значению за ограниченное число шагов. Поскольку вычисление и стоит флопов, для вычисления одного собственного значения достаточно , а для вычисления всех собственных значений матрицы требуется флопов. Для собственных векторов матрицы мы легко можем получить явные выражения.
Лемма 2. Если - собственное значение матрицы , то соответствующий вектор равен . Поскольку матрица диагональная, для вычисления такого вектора достаточно флопов.
Доказательство.
поскольку , что и требовалось.
Для вычисления по этой простой формуле всех собственных векторов требуется флопов. К сожалению, формула не обеспечивает численной устойчивости, так как для двух очень близких значений может давать неортогональные приближенные собственные векторы . Потребовалось целое десятилетие для того, чтобы найти устойчивую альтернативу исходному описанию алгоритма. Снова детали будут обсуждаться позднее в данном разделе. Алгоритм является рекурсивным. 1.2.1 Дефляция[править] До сих пор полагалось, что все различны и все отличны от нуля. Если это не так, вековое уравнение имеет вертикальных асимптот, где , а потому корней. Однако оказывается, что остальные собственных значений могут быть определены без каких-либо усилий: если или , то легко показать, что является собственным значением и для матрицы . В такой ситуации мы говорим о дефляции. На практике выбирается некоторое пороговое значение и дефляция для числа регистрируется, если в смысле этого порога достаточно близко к либо достаточно мало. Основной выигрыш от использования дефляции состоит не в том, что убыстряется решение векового уравнения - этот этап в любом случае стоит лишь операций. Выигрыш заключается в ускорении матричного умножения на последнем шаге алгоритма. Действительно, если , то соответствующий собственный вектор есть i-й столбец единичной матрицы. Это означает, что является i-м столбцом в матрице , поэтому при формировании матрицы посредством левого умножения на вычисление i-го столбца не требует никаких затрат. Аналогичное упрощение имеет место в случае . При дефляции многих собственных значений устраняется большая часть работы, связанной с матричным умножением. 1.3 Вычислительное ядро алгоритма[править] Вычислительным ядром последовательной схемы решения является вычисление матрицы собственных векторов путём умножения матрицы на матрицу Данная операция имеет сложность о чём и говорится в разделе#Последовательная сложность алгоритма . Ей предшествует вычисление собственных значений и векторов матрицы 1.4 Макроструктура алгоритма[править] В разделе #Информационный граф описана структура алгоритма, в которой есть блок умножения матриц для вычисления собственных векторов, являющийся вычислительным ядром алгоритма. В соответствующем разделе (#Вычислительное ядро алгоритма) мы упоминали о том, что данному блоку предшествует вычисление собственных значений, которое производится методом Ньютона. 1.4.1 Решение векового уравнения[править] Подробно стоит поговорить о решении векового уравнения, которое является одной из основных частей алгоритма. Предположим, что некоторое , хотя и мало, все же недостаточно мало для того, чтобы была зарегистрирована дефляция. В этом случае применение метода Ньютона к решению векового уравнения встречается с затруднениями. Вспомним, что пересчет приближённого решения уравнения в методе Ньютона основан на следующих положениях: 1. Вблизи точки функция аппроксимируется линейной функцией ; график есть прямая линия, касающаяся графика функции при . 2. В качестве берётся нуль этого линейного приближения, т.е. . Функция, показанная на рис.1, не доставляет видимых трудностей методу Ньютона, поскольку вблизи каждого своего нуля достаточно хорошо аппроксимируется линейными функциями. Однако рассмотрим график функции на рис. 2. Она получена из функции на рис. 1 заменой значения .5 для на .001. Это новое значение недостаточно мало для того, чтобы вызвать дефляцию. График функции в левой части рис.2 визуально не отличим от её вертикальных и горизонтальных асимптот, поэтому в правой части укрупненно воспроизведён фрагмент графика, прилегающий к вертикальной асимптоте . Видно, что график слишком быстро "выполняет поворот" и для большей части значений почти горизонтален. Поэтому, применяя метод Ньютона почти к любому начальному приближению , мы получаем линейное приближение с почти горизонтальным графиком и малым положительным угловым коэффициентом. В результате является отрицательным числом, огромным по абсолютной величине, которое совершенно бесполезно в качестве приближения к истинному корню.
Рис. 2. График функции Чтобы найти выход из этого положения, можно модифицировать метод Ньютона следующим образом: раз нельзя хорошо приблизить линейной функцией , попробуем взять в качестве приближения какую-нибудь другую простую функцию . Нет ничего особого именно в прямых линиях: для метода Ньютона вместо можно взять любое приближение , значения и нули которого легко вычисляются. Функция имеет полюсы в точках и , которые определяют её поведение в соответствующих окрестностях. Поэтому при поиске корня в интервале естественно выбрать функцию , также имеющую эти полюсы, т.е. функцию вида
Константы и обеспечивающие, что есть приближение к , можно выбрать несколькими способами. Отметим, что если и уже известны, то уравнение легко решается относительно , поскольку сводится к эквивалентному квадратному уравнению
Пусть - приближённое значение корня. определим и так, чтобы
для в окрестности . Заметим, что
.
Если , то есть сумма положительных слагаемых, а - сумма отрицательных. Поэтому и , и могут быть вычислены с высокой точностью; однако при их сложении вполне вероятно взаимное уничтожение верных разрядов и потеря относительной точности в сумме. Возьмем числа и , такие, что функция
удовлетворяет условиям и (1)
Это означает, что гипербола, являющаяся графиком функции , касается графика функции при . Два условия в (1) - это обычные условия метода Ньютона, за исключением того, что вместо прямой в качестве приближения используется гипербола. Легко проверить, что и Подобным же образом выбираем и так, чтобы функция
Удовлетворяла условиям
и (2)
Наконец, полагаем
.
1.5 Схема реализации последовательного алгоритма[править] Вычисление собственных значений и собственных векторов симметричной трехдиагональной матрицы посредством стратегии "разделяй и властвуй":
proc dc_eig по входной матрице вычисляются выходные матрицы и , такие, что Если - матрица размера x
1. Присвоить выходным параметрам значения
Иначе
1. Представить в виде 2. call dc_eig 3. call dc_eig 4. Построить по 5. Найти матрицу собственных значений 6. Найти матрицу собственных векторов для матрицы 7. Построить матрицу собственных векторов для матрицы : 8. Присвоить выходным параметрам значения и
endif 1.6 Последовательная сложность алгоритма[править] Проанализируем сложность алгоритма. Пусть - число флопов при обработке матрицы размера процедурой dc_eig. Тогда
два рекурсивных обращения к dc_eig вычисление собственных значений матрицы вычисление собственных векторов матрицы вычисление матрицы
Если и рассматриваются как плотные матрицы и используется стандартный алгоритм матричного умножения, то константа в последней строке равна 1. Таким образом, именно это умножение составляет наиболее трудоёмкую часть алгоритма в целом. Игнорируя члены порядка , получаем . Решая это разностное уравнение, находим На практике константа обычно гораздо меньше 1, потому что матрица весьма разрежена вследствие явления, называемого дефляцией. 1.7 Информационный граф[править] В данном разделе представлен информационный граф алгоритма: на рисунке 2 изображена структура всего алгоритма в целом, в то время как на рисунке 3 детально описана одна из ячеек структуры.
Рис. 3. Дерево алгоритма "Разделяй и властвуй"
Рис. 4. Детальное описание одного блока алгоритма 1.8 Ресурс параллелизма алгоритма[править] Рассмотрев информационный граф алгоритма, можно заметить, что в структуре имеется два параллельных блока - вызов рекурсивных функций dc_eig для вычисления собственных значений и векторов матриц и - это единственная часть алгоритма, в которой мы прибегаем к параллелизму. В случае реализации без параллелизма, функции dc_eig отрабатывают последовательно - сначала со входными параметрами , затем - . 1.9 Входные и выходные данные алгоритма[править] На вход алгоритму подаётся трёхдиагональная матрица, описанная в разделе #Математическое описание алгоритма На выходе мы получаем собственные значения и собственные вектора исходной матрицы 1.10 Свойства алгоритма[править] Описанный в данной статье алгоритм является самым быстрым алгоритмом, среди существующих: QR / Бисекция и обратная итерация / Разделяй и властвуй. Последовательная сложность алгоритма: Параллельная сложность алгоритма (в силу использования дефляции): или в самых редких случаях (см. [3] стр. 8). Особенности алгоритма: 1. Использование дефляции (#Дефляция) 2. Использование адаптированного метода Ньютона (#Макроструктура алгоритма) 2 ЧАСТЬ. Программная реализация алгоритма[править] Вторая часть описания алгоритмов в рамках AlgoWiki рассматривает все составные части процесса их реализации. Рассматривается как последовательная реализация алгоритма, так и параллельная. Описывается взаимосвязь свойств программ, реализующих алгоритм, и особенностей архитектуры компьютера, на которой они выполняются. Исследуется работа с памятью, локальность данных и вычислений, описывается масштабируемость и эффективность параллельных программ, производительность компьютеров, достигаемая на данной программе. Обсуждаются особенности реализации для разных классов архитектур компьютеров, приводятся ссылки на реализации в существующих библиотеках.
2.1 Масштабируемость алгоритма и его реализации[править] TBD
2.2 Существующие реализации алгоритма[править] Большого количества реализаций данного алгоритма найдено не было: единственный пример описан на страницах 7, 19 и 20 в этой статье, которая указана в списке литературы.
