Метод "Разделяй и властвуй" Завольсков/Землянский: различия между версиями
| [непроверенная версия] | [непроверенная версия] |
| Строка 85: | Строка 85: | ||
т.е. собственные значения матрицы <math>T</math> есть корни так называемого векового уравнения <math>f(\lambda) = 0</math>. Если все числа <math>d_{i}</math> различны и все <math>u_{i} <> 0</math> (случай общего положения), то <math>f(\lambda)</math> имеет график типа показанного на рис.<math> 1 </math>(где <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> имеет график типа показанного на рис.<math> 1 </math>(где <math>n = 4</math> и <math>\rho > 0</math>). | ||
[[File:Graphics.PNG|1000px|Рис. 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:Graphics.PNG|1000px|thumb|Рис. 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>]] | ||
Можно видеть, что прямая <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> флопов. | Можно видеть, что прямая <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> флопов. | ||
Версия от 15:28, 14 октября 2016
Авторы статьи: Завольсков Владислав, Землянский Роман, 614 группа
ЧАСТЬ. Свойства и структура алгоритмов
Общее описание алгоритма
Метод разделяй и властвуй вычисления собственных значений и векторов трёхдиагональной матрицы - это наиболее быстрый из существующих методов, если нужны все собственные значения и собственные векторы трехдиагональной матриц, начиная с порядка n, примерно равного 26. (Точное значение этого порогового порядка зависит от компьютера.) Его численно устойчивая реализация весьма не тривиальна. В самом деле, хотя впервые метод был предложен еще в 1981 г., "правильный" способ его реализации был найден лишь в 1992 г. Этот способ воплощен LAPACK-программами ssyevd (для плотных матриц) и sstevd (для трехдиагональных матриц). В них стратегия "разделяй-и-влавствуй" используется для матриц порядка, большего чем 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} + \begin{bmatrix} &&&&& \\ &&b_{m} & b_{m} \\ &&b_{m} & b_{m} \\ &&&&& \\ \end{bmatrix} = }[/math]
- [math]\displaystyle{ = \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} = \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{ v = \begin{bmatrix} 0 , \ldots , 0 , 1 , 1 , 0 \ldots , 0 \end{bmatrix}^T }[/math], получим матрицу, состоящую из последнего столбца матрицы [math]\displaystyle{ Q_{1}^{T} }[/math] и первого столбца матрицы [math]\displaystyle{ Q_{2}^{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} = 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] имеет график типа показанного на рис.[math]\displaystyle{ 1 }[/math](где [math]\displaystyle{ n = 4 }[/math] и [math]\displaystyle{ \rho \gt 0 }[/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 }[/math]
[math]\displaystyle{ =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 }[/math]
поскольку [math]\displaystyle{ \rho u^{T}(D - \alpha I)^{-1}u + 1 = f(\alpha) = 0 }[/math]
[math]\displaystyle{ =\alpha [(D - \alpha I)^{-1}u] }[/math], что и требовалось.
Для вычисления по этой простой формуле всех [math]\displaystyle{ n }[/math] собственных векторов требуется [math]\displaystyle{ O(n^{2}) }[/math] флопов. К сожалению, формула не обеспечивает численной устойчивости, так как для двух очень близких значений [math]\displaystyle{ \alpha_{i} }[/math] может давать неортогональные приближенные собственные векторы [math]\displaystyle{ u_{i} }[/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 x 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 x n }[/math] процедурой dc_eig. Тогда
[math]\displaystyle{ t(n) = 2t(n/2) }[/math] два рекурсивных обращения к dc_eig<(math>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 = c\frac{4}{3}n^{3} }[/math]
На практике константа [math]\displaystyle{ c }[/math] обычно гораздо меньше 1, потому что матрица [math]\displaystyle{ Q^{'} }[/math] весьма разрежена вследствие явления, называемого дефляцией.
Информационный граф
TBD
Ресурс параллелизма алгоритма
TBD
Входные и выходные данные алгоритма
В данном разделе необходимо описать объем, структуру, особенности и свойства входных и выходных данных алгоритма: векторы, матрицы, скаляры, множества, плотные или разреженные структуры данных, их объем. Полезны предположения относительно диапазона значений или структуры, например, диагональное преобладание в структуре входных матриц, соотношение между размером матриц по отдельным размерностям, большое число матриц очень малой размерности, близость каких-то значений к машинному нулю, характер разреженности матриц и другие.
Свойства алгоритма
Какая-нибудь шляпа
ЧАСТЬ. Программная реализация алгоритма
Вторая часть описания алгоритмов в рамках AlgoWiki рассматривает все составные части процесса их реализации. Рассматривается как последовательная реализация алгоритма, так и параллельная. Описывается взаимосвязь свойств программ, реализующих алгоритм, и особенностей архитектуры компьютера, на которой они выполняются. Исследуется работа с памятью, локальность данных и вычислений, описывается масштабируемость и эффективность параллельных программ, производительность компьютеров, достигаемая на данной программе. Обсуждаются особенности реализации для разных классов архитектур компьютеров, приводятся ссылки на реализации в существующих библиотеках.
Масштабируемость алгоритма и его реализации
TBD
Существующие реализации алгоритма
TBD
Литература
[1] Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. - СПб.: БХВ-Петербург, 2002. - 608 с.
[2] Воеводин В.В., Воеводин Вад.В. Спасительная локальность суперкомпьютеров //Открытые системы. - 2013. - № 9. - С. 12-15.
[3] Воеводин Вад.В., Швец П. Метод покрытий для оценки локальности использования данных в программах // Вестник УГАТУ. — 2014. — Т. 18, № 1(62). — С. 224–229.
[4] Антонов А.С., Теплов А.М. О практической сложности понятия масштабируемости параллельных программ// Высокопроизводительные параллельные вычисления на кластерных системах (HPC 2014): Материалы XIV Международной конференции -Пермь: Издательство ПНИПУ, 2014. С. 20-27.
[5] Никитенко Д.А. Комплексный анализ производительности суперкомпьютерных систем, основанный на данных системного мониторинга // Вычислительные методы и программирование. 2014. 15. 85–97.