|
|
| Строка 2: |
Строка 2: |
| | | | |
| | = ЧАСТЬ. Свойства и структура алгоритмов = | | = ЧАСТЬ. Свойства и структура алгоритмов = |
| − | == Общее описание алгоритма ==
| |
| − | '''Метод разделяй и властвуй вычисления собственных значений и векторов трёхдиагональной матрицы''' - это наиболее быстрый из существующих методов, если нужны все собственные значения и собственные векторы трехдиагональной матриц, начиная с порядка n, примерно равного 26. (Точное значение этого порогового порядка зависит от компьютера.) Его численно устойчивая реализация весьма не тривиальна. В самом деле, хотя впервые метод был предложен еще в 1981 г., "правильный" способ его реализации был найден лишь в 1992 г. Этот способ воплощен LAPACK-программами ssyevd (для плотных матриц) и sstevd (для трехдиагональных матриц). В них стратегия "разделяй-и-влавствуй" используется для матриц порядка, большего чем 25. Для матриц меньшего порядка (или если нужны только собственные значения) происходит автоматический переход к QR-итерации.
| |
| − |
| |
| − | == Математическое описание алгоритма ==
| |
| − | Пусть
| |
| − |
| |
| − | :<math>
| |
| − | 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>
| |
| − | + \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}
| |
| − |
| |
| − | = \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T}
| |
| − | </math>
| |
| − |
| |
| − | Предположим, что нам известны спектральные разложения матриц <math>T_{1}</math> и <math> T_{2} </math>:
| |
| − | <math> T_{i} = Q_{i} \Lambda_{i} Q_{i}^{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}vv^{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>v = \begin{bmatrix} 0 , \ldots , 0 , 1 , 1 , 0 \ldots , 0 \end{bmatrix}^T</math>, получим матрицу, состоящую из последнего столбца матрицы <math> Q_{1}^{T}</math> и первого столбца матрицы <math> Q_{2}^{T}</math>.
| |
| − |
| |
| − | Следовательно, <math>T</math> имеет те же собственные значения, что и пдобная ей матрица <math>D + \rho uu^{T}</math>,
| |
| − | где <math>D = \begin{bmatrix} L_{1} & 0 \\ 0 & L_{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} <= \ldots <=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} = f(\lambda)</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|thumb|center|800px|Рис. 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>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</math>
| |
| − |
| |
| − | <math>=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</math>
| |
| − |
| |
| − | поскольку <math> \rho u^{T}(D - \alpha I)^{-1}u + 1 = f(\alpha) = 0 </math>
| |
| − |
| |
| − | <math>=\alpha [(D - \alpha I)^{-1}u]</math>, что и требовалось.
| |
| − |
| |
| − |
| |
| − | Для вычисления по этой простой формуле всех <math>n</math> собственных векторов требуется <math>O(n^{2})</math> флопов. К сожалению, формула не обеспечивает численной устойчивости, так как для двух очень близких значений <math>\alpha_{i}</math> может давать неортогональные приближенные собственные векторы <math>u_{i}</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>
| |
| − |
| |
| − | Если <math>T</math> - матрица размера <math>1 x 1</math>
| |
| − |
| |
| − | 1. Присвоить выходным параметрам значения <math> Q = 1, \Lambda = T</math>
| |
| − |
| |
| − | Иначе
| |
| − |
| |
| − | 1. Представить <math>T</math> в виде <math> T = \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T} </math>
| |
| − |
| |
| − | 2. ''call dc_eig''<math>(T_{1},Q_{1},\Lambda_{1})</math>
| |
| − |
| |
| − | 3. ''call dc_eig''<math>(T_{2},Q_{2},\Lambda_{2})</math>
| |
| − |
| |
| − | 4. Построить <math>D+\rho uu^{T}</math> по <math> Lambda_{1},Lambda_{2}, Q_{1}, Q_{2}</math>
| |
| − |
| |
| − | 5. Найти матрицу собственных значений <math>\Lambda</math>
| |
| − |
| |
| − | 6. Найти матрицу собственных векторов <math>Q^{'}</math> для матрицы <math>D+\rho uu^{T}</math>
| |
| − |
| |
| − | 7. Построить матрицу собственных векторов <math>Q</math> для матрицы <math>T</math> :
| |
| − |
| |
| − | <math> Q = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}* Q^{'} </math>
| |
| − |
| |
| − | 8. Присвоить выходным параметрам значения <math>Q</math> и <math>\Lambda</math>
| |
| − |
| |
| − | ''endif''
| |
| − |
| |
| − | == Последовательная сложность алгоритма ==
| |
| − | Проанализируем сложность алгоритма. Пусть <math>t(n)</math> - число флопов при обработке матрицы размера <math> n x 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 = c\frac{4}{3}n^{3} </math>
| |
| − |
| |
| − | На практике константа <math>c</math> обычно гораздо меньше 1, потому что матрица <math>Q^{'}</math> весьма разрежена вследствие явления, называемого дефляцией.
| |
| − |
| |
| − | == Информационный граф ==
| |
| | TBD | | TBD |
| − |
| |
| − | == Ресурс параллелизма алгоритма ==
| |
| − | TBD
| |
| − |
| |
| − | == Входные и выходные данные алгоритма ==
| |
| − | В данном разделе необходимо описать объем, структуру, особенности и свойства входных и выходных данных алгоритма: векторы, матрицы, скаляры, множества, плотные или разреженные структуры данных, их объем. Полезны предположения относительно диапазона значений или структуры, например, диагональное преобладание в структуре входных матриц, соотношение между размером матриц по отдельным размерностям, большое число матриц очень малой размерности, близость каких-то значений к машинному нулю, характер разреженности матриц и другие.
| |
| − |
| |
| − | == Свойства алгоритма ==
| |
| − | Какая-нибудь шляпа
| |
| − |
| |
| − | = ЧАСТЬ. Программная реализация алгоритма =
| |
| − | Вторая часть описания алгоритмов в рамках AlgoWiki рассматривает все составные части процесса их реализации. Рассматривается как последовательная реализация алгоритма, так и параллельная. Описывается взаимосвязь свойств программ, реализующих алгоритм, и особенностей архитектуры компьютера, на которой они выполняются. Исследуется работа с памятью, локальность данных и вычислений, описывается масштабируемость и эффективность параллельных программ, производительность компьютеров, достигаемая на данной программе. Обсуждаются особенности реализации для разных классов архитектур компьютеров, приводятся ссылки на реализации в существующих библиотеках.
| |
| − |
| |
| − |
| |
| − | == Масштабируемость алгоритма и его реализации ==
| |
| − | TBD
| |
| − |
| |
| − |
| |
| − | == Существующие реализации алгоритма ==
| |
| − | TBD
| |
| − |
| |
| | = Литература = | | = Литература = |
| | [1] Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. - СПб.: БХВ-Петербург, 2002. - 608 с. | | [1] Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. - СПб.: БХВ-Петербург, 2002. - 608 с. |
