Приложение 1: различия между версиями
[непроверенная версия] | [непроверенная версия] |
ASA (обсуждение | вклад) (Полностью удалено содержимое страницы) |
ASA (обсуждение | вклад) |
||
Строка 1: | Строка 1: | ||
+ | = Разложение Холецкого (метод квадратного корня) = | ||
+ | == Свойства и структура алгоритма == | ||
+ | |||
+ | === Общее описание алгоритма === | ||
+ | |||
+ | '''Разложение Холецкого''' впервые предложено французским офицером и математиком Андре-Луи Холецким в конце Первой Мировой войны, незадолго до его гибели в бою в августе 1918 г. Идея этого разложения была опубликована в 1924 г. его сослуживцем<ref>Commandant Benoit, Note sur une méthode de résolution des équations normales provenant de l'application de la méthode des moindres carrés à un système d'équations linéaires en nombre inférieur à celui des inconnues (Procédé du Commandant Cholesky), Bulletin Géodésique 2 (1924), 67-77.</ref>. Потом оно было использовано поляком Т. Банашевичем<ref>Banachiewicz T. Principles d'une nouvelle technique de la méthode des moindres carrês. Bull. Intern. Acad. Polon. Sci. A., 1938, 134-135.</ref><ref>Banachiewicz T. Méthode de résoltution numérique des équations linéaires, du calcul des déterminants et des inverses et de réduction des formes quardatiques. Bull. Intern. Acad. Polon. Sci. A., 1938, 393-401.</ref> в 1938 г. В советской математической литературе называется также методом квадратного корня<ref>Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977.</ref><ref>Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. М.: Наука, 1984.</ref><ref>Фаддеев Д.К., Фаддева В.Н. Вычислительные основы линейной алгебры. М.-Л.: Физматгиз, 1963.</ref>; название связано с характерными операциями, отсутствующими в родственном разложении Гаусса. | ||
+ | |||
+ | Первоначально разложение Холецкого использовалось исключительно для плотных симметричных положительно определенных матриц. В настоящее время его использование гораздо шире. Оно может быть применено также, например, к эрмитовым матрицам. Для повышения производительности вычислений часто применяется блочная версия разложения. | ||
+ | |||
+ | Для разреженных матриц разложение Холецкого также широко применяется в качестве основного этапа прямого метода решения линейных систем. В этом случае используют специальные упорядочивания для уменьшения ширины профиля исключения, а следовательно и уменьшения количества арифметических операций. Другие упорядочивания используются для выделения независимых блоков вычислений при работе на системах с параллельной организацией. | ||
+ | |||
+ | Варианты разложения Холецкого нашли успешные применения и в итерационных методах для построения переобусловливателей разреженных симметричных положительно определенных матриц. В неполном треугольном разложении («по позициям») элементы переобусловливателя вычисляются только в заранее заданных позициях, например, в позициях ненулевых элементов исходной матрицы (так называемое разложение IC0). Для получения же более точного разложения применяется приближение, в котором фильтрация малых элементов производится «по значениям». В зависимости от используемого порога фильтрации можно получить более точное, но и более заполненное разложение. Существует и алгоритм разложения второго порядка точности<ref>Kaporin I.E. High quality preconditioning of a general symmetric positive definite matrix based on its UTU + UTR + RTU-decomposition. Numer. Lin. Algebra Appl. (1998) Vol. 5, No. 6, 483-509.</ref>. В нём при таком же заполнении множителей разложения удается улучшить точность. Для такого разложения в параллельном режиме используется специальный вариант аддитивного переобуславливания на основе разложения второго порядка<ref>Капорин И.Е., Коньшин И.Н. Параллельное решение симметричных положительно-определенных систем на основе перекрывающегося разбиения на блоки. Ж. вычисл. матем. и матем. физ., 2001, Т, 41, N. 4, C. 515–528.</ref>. | ||
+ | |||
+ | На этой странице представлено исходное разложение Холецкого с новых позиций нашего суперкомпьютерного века. Приведено описание конкретной версии разложения Холецкого — для плотных вещественных симметричных положительно определённых матриц, но структура для ряда других версий, например, для комплексного случая, почти такая же, различия состоят в замене большинства вещественных операций на комплексные. Список остальных основных вариантов разложения можно посмотреть на странице [[Метод Холецкого (нахождение симметричного треугольного разложения)]]. | ||
+ | |||
+ | Используется для разложения положительно определённых эрмитовых (''в вещественном случае - симметрических'') матриц в виде <math>A = L L^*</math> (<math>L</math> — [http://lineal.guru.ru/lineal3/texts.php?PHPSESSID=e8862lf95ipuht6hb7bd6vl1m6&&rd=3368 нижняя треугольная матрица]) или в виде <math>A = U^* U</math> (<math>U</math> — [http://lineal.guru.ru/lineal3/texts.php?PHPSESSID=e8862lf95ipuht6hb7bd6vl1m6&&rd=3368 верхняя треугольная матрица] ; эти разложения совершенно эквивалентны друг другу по вычислениям и различаются только способом представления данных). Он заключается в реализации формул для элементов матрицы <math>L</math>, получающихся из вышеприведённого равенства единственным образом. Получило широкое распространение благодаря следующим особенностям. | ||
+ | |||
+ | ==== Симметричность и положительная определённость матрицы ==== | ||
+ | |||
+ | Симметричность матрицы позволяет хранить и вычислять только чуть больше половины её элементов, что почти вдвое экономит как необходимые для вычислений объёмы памяти, так и количество операций в сравнении с, например, разложением по методу Гаусса. При этом альтернативное (без вычисления квадратных корней) LU-разложение, использующее симметрию матрицы, всё же несколько быстрее метода Холецкого (не использует извлечение квадратных корней), но требует хранения всей матрицы. | ||
+ | Благодаря тому, что разлагаемая матрица не только симметрична, но и положительно определена, её LU-разложения, в том числе и разложение методом Холецкого, имеют наименьшее ''[[Глоссарий#Эквивалентное возмущение|эквивалентное возмущение]]'' из всех известных разложений матриц. | ||
+ | |||
+ | ==== Режим накопления ==== | ||
+ | |||
+ | Алгоритм позволяет использовать так называемый ''режим накопления'', обусловленный тем, что значительную часть вычислений составляют ''вычисления скалярных произведений''. | ||
+ | |||
+ | === Математическое описание алгоритма === | ||
+ | |||
+ | Исходные данные: положительно определённая симметрическая матрица <math>A</math> (элементы <math>a_{ij}</math>). | ||
+ | |||
+ | Вычисляемые данные: нижняя треугольная матрица <math>L</math> (элементы <math>l_{ij}</math>). | ||
+ | |||
+ | Формулы метода: | ||
+ | :<math> | ||
+ | \begin{align} | ||
+ | l_{11} & = \sqrt{a_{11}}, \\ | ||
+ | l_{j1} & = \frac{a_{j1}}{l_{11}}, \quad j \in [2, n], \\ | ||
+ | l_{ii} & = \sqrt{a_{ii} - \sum_{p = 1}^{i - 1} l_{ip}^2}, \quad i \in [2, n], \\ | ||
+ | l_{ji} & = \left (a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp} \right ) / l_{ii}, \quad i \in [2, n - 1], j \in [i + 1, n]. | ||
+ | \end{align} | ||
+ | </math> | ||
+ | |||
+ | Существует также блочная версия метода, однако в данном описании разобран только точечный метод. | ||
+ | |||
+ | В ряде реализаций деление на диагональный элемент выполняется в два этапа: вычисление <math>\frac{1}{l_{ii}}</math> и затем умножение на него всех (видоизменённых) <math>a_{ji}</math> . Здесь мы этот вариант алгоритма не рассматриваем. Заметим только, что он имеет худшие параллельные характеристики, чем представленный. | ||
+ | |||
+ | === Вычислительное ядро алгоритма === | ||
+ | |||
+ | Вычислительное ядро последовательной версии метода Холецкого можно составить из множественных (всего их <math>\frac{n (n - 1)}{2}</math>) вычислений скалярных произведений строк матрицы: | ||
+ | |||
+ | :<math>\sum_{p = 1}^{i - 1} l_{ip} l_{jp}</math> | ||
+ | |||
+ | в режиме накопления или без него, в зависимости от требований задачи. Во многих последовательных реализациях упомянутый способ представления не используется. Дело в том, что в них вычисление сумм типа | ||
+ | |||
+ | :<math>a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp}</math><nowiki/> | ||
+ | |||
+ | в которых и встречаются скалярные произведения, ведутся не в порядке «вычислили скалярное произведение, а потом вычли его из элемента», а путём вычитания из элемента покомпонентных произведений, являющихся частями скалярных произведений. Поэтому следует считать вычислительным ядром метода не вычисления скалярных произведений, а вычисления выражений | ||
+ | |||
+ | :<math>a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp}</math> | ||
+ | |||
+ | в ''режиме накопления'' или без него. | ||
+ | |||
+ | Тем не менее, в популярных зарубежных реализациях точечного метода Холецкого, в частности, в библиотеках LINPACK и LAPACK, основанных на BLAS, используются именно вычисления скалярных произведений в виде вызова соответствующих подпрограмм BLAS (конкретно — функции SDOT). На последовательном уровне это влечёт за собой дополнительную операцию суммирования на каждый из <math>\frac{n (n + 1)}{2}</math> вычисляемый элемент матрицы <math>L</math> и некоторое замедление работы программы (о других следствиях рассказано ниже в разделе «[[#Существующие реализации алгоритма|Существующие реализации алгоритма]]»). Поэтому в данных вариантах ядром метода Холецкого будут вычисления этих скалярных произведений. | ||
+ | |||
+ | === Макроструктура алгоритма === | ||
+ | |||
+ | Как записано и в [[#Вычислительное ядро алгоритма|описании ядра алгоритма]], основную часть метода составляют множественные (всего <math>\frac{n (n - 1)}{2}</math>) вычисления сумм | ||
+ | |||
+ | :<math>a_{ji}-\sum_{p=1}^{i-1}l_{ip} l_{jp}</math> | ||
+ | |||
+ | в режиме накопления или без него. | ||
+ | |||
+ | === Схема реализации последовательного алгоритма === | ||
+ | |||
+ | Последовательность исполнения метода следующая: | ||
+ | |||
+ | 1. <math>l_{11}= \sqrt{a_{11}}</math> | ||
+ | |||
+ | 2. <math>l_{j1}= \frac{a_{j1}}{l_{11}}</math> (при <math>j</math> от <math>2</math> до <math>n</math>). | ||
+ | |||
+ | Далее для всех <math>i</math> от <math>2</math> до <math>n</math> по нарастанию выполняются | ||
+ | |||
+ | 3. <math>l_{ii} = \sqrt{a_{ii} - \sum_{p = 1}^{i - 1} l_{ip}^2}</math> и | ||
+ | |||
+ | 4. (кроме <math>i = n</math>): <nowiki/><math>l_{ji} = \left (a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp} \right ) / l_{ii}</math> (для всех <math>j</math> от <math>i + 1</math> до <math>n</math>). | ||
+ | |||
+ | После этого (если <math>i < n</math>) происходит переход к шагу 3 с бо́льшим <math>i</math>. | ||
+ | |||
+ | Особо отметим, что вычисления сумм вида <math>a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp}</math> в обеих формулах производят в режиме накопления вычитанием из <math>a_{ji}</math> произведений <math>l_{ip} l_{jp}</math> для <math>p</math> от <math>1</math> до <math>i - 1</math>, c нарастанием <math>p</math>. | ||
+ | |||
+ | === Последовательная сложность алгоритма === | ||
+ | |||
+ | Для разложения матрицы порядка n методом Холецкого в последовательном (наиболее быстром) варианте требуется: | ||
+ | |||
+ | * <math>n</math> вычислений квадратного корня, | ||
+ | * <math>\frac{n(n-1)}{2}</math> делений, | ||
+ | * <math>\frac{n^3-n}{6}</math> сложений (вычитаний), | ||
+ | * <math>\frac{n^3-n}{6}</math> умножений. | ||
+ | |||
+ | Умножения и сложения (вычитания) составляют ''основную часть алгоритма''. | ||
+ | |||
+ | При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности (или использования функции вроде DPROD в Фортране), что ещё больше увеличивает долю умножений и сложений/вычитаний во времени, требуемом для выполнения метода Холецкого. | ||
+ | |||
+ | При классификации по последовательной сложности, таким образом, метод Холецкого относится к алгоритмам ''с кубической сложностью''. | ||
+ | |||
+ | === Информационный граф === | ||
+ | |||
+ | Опишем [[глоссарий#Граф алгоритма|граф алгоритма]]<ref>Воеводин В.В. Математические основы параллельных вычислений// М.: Изд. Моск. ун-та, 1991. 345 с.</ref><ref>Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. – СПб.: БХВ - Петербург, 2002. – 608 с.</ref><ref>Фролов А.В.. Принципы построения и описание языка Сигма. Препринт ОВМ АН N 236. М.: ОВМ АН СССР, 1989.</ref> как аналитически, так и в виде рисунка. | ||
+ | |||
+ | Граф алгоритма состоит из трёх групп вершин, расположенных в целочисленных узлах трёх областей разной размерности. | ||
+ | |||
+ | '''Первая''' группа вершин расположена в одномерной области, соответствующая ей операция вычисляет функцию SQRT. | ||
+ | Единственная координата каждой из вершин <math>i</math> меняется в диапазоне от <math>1</math> до <math>n</math>, принимая все целочисленные значения. | ||
+ | |||
+ | Аргумент этой функции | ||
+ | |||
+ | * при <math>i = 1</math> — элемент ''входных данных'', а именно <math>a_{11}</math>; | ||
+ | * при <math>i > 1</math> — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами <math>i - 1</math>, <math>i</math>, <math>i - 1</math>. | ||
+ | Результат срабатывания операции является ''выходным данным'' <math>l_{ii}</math>. | ||
+ | |||
+ | '''Вторая''' группа вершин расположена в двумерной области, соответствующая ей операция <math>a / b</math>. | ||
+ | Естественно введённые координаты области таковы: | ||
+ | * <math>i</math> — меняется в диапазоне от <math>1</math> до <math>n-1</math>, принимая все целочисленные значения; | ||
+ | * <math>j</math> — меняется в диапазоне от <math>i+1</math> до <math>n</math>, принимая все целочисленные значения. | ||
+ | |||
+ | Аргументы операции следующие: | ||
+ | *<math>a</math>: | ||
+ | ** при <math>i = 1</math> — элементы ''входных данных'', а именно <math>a_{j1}</math>; | ||
+ | ** при <math>i > 1</math> — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами <math>i - 1, j, i - 1</math>; | ||
+ | * <math>b</math> — результат срабатывания операции, соответствующей вершине из первой группы, с координатой <math>i</math>. | ||
+ | |||
+ | Результат срабатывания операции является ''выходным данным'' <math>l_{ji}</math>. | ||
+ | |||
+ | '''Третья''' группа вершин расположена в трёхмерной области, соответствующая ей операция <math>a - b * c</math>. | ||
+ | Естественно введённые координаты области таковы: | ||
+ | * <math>i</math> — меняется в диапазоне от <math>2</math> до <math>n</math>, принимая все целочисленные значения; | ||
+ | * <math>j</math> — меняется в диапазоне от <math>i</math> до <math>n</math>, принимая все целочисленные значения; | ||
+ | * <math>p</math> — меняется в диапазоне от <math>1</math> до <math>i - 1</math>, принимая все целочисленные значения. | ||
+ | |||
+ | Аргументы операции следующие: | ||
+ | *<math>a</math>: | ||
+ | ** при <math>p = 1</math> элемент входных данных <math>a_{ji}</math>; | ||
+ | ** при <math>p > 1</math> — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами <math>i, j, p - 1</math>; | ||
+ | *<math>b</math> — результат срабатывания операции, соответствующей вершине из второй группы, с координатами <math>p, i</math>; | ||
+ | *<math>c</math> — результат срабатывания операции, соответствующей вершине из второй группы, с координатами <math>p, j</math>; | ||
+ | |||
+ | Результат срабатывания операции является ''промежуточным данным'' алгоритма. | ||
+ | |||
+ | Описанный граф можно посмотреть на рис.1 и рис.2, выполненных для случая <math>n = 4</math>. Здесь вершины первой группы обозначены жёлтым цветом и буквосочетанием sq, вершины второй — зелёным цветом и знаком деления, третьей — красным цветом и буквой f. Вершины, соответствующие операциям, производящим выходные данные алгоритма, выполнены более крупно. Дублирующие друг друга дуги даны как одна. На рис.1 показан граф алгоритма согласно классическому определению , на рис.2 к графу алгоритма добавлены вершины , соответствующие входным (обозначены синим цветом) и выходным (обозначены розовым цветом) данным. | ||
+ | |||
+ | [[file:Cholesky full.png|thumb|center|1400px|Рисунок 1. Граф алгоритма без отображения входных и выходных данных. SQ - вычисление квадратного корня, F - операция a-bc, Div - деление.]] | ||
+ | [[file:Cholesky full_in_out.png|thumb|center|1400px|Рисунок 2. Граф алгоритма с отображением входных и выходных данных. SQ - вычисление квадратного корня, F - операция a-bc, Div - деление, In - входные данные, Out - результаты.]] | ||
+ | |||
+ | === Ресурс параллелизма алгоритма === | ||
+ | |||
+ | Для разложения матрицы порядка <math>n</math> методом Холецкого в параллельном варианте требуется последовательно выполнить следующие ярусы: | ||
+ | * <math>n</math> ярусов с вычислением квадратного корня (единичные вычисления в каждом из ярусов), | ||
+ | * <math>n - 1</math> ярус делений (в каждом из ярусов линейное количество делений, в зависимости от яруса — от <math>1</math> до <math>n - 1</math>), | ||
+ | * по <math>n - 1</math> ярусов умножений и сложений/вычитаний (в каждом из ярусов — квадратичное количество операций, от <math>1</math> до <math>\frac{n^2 - n}{2}</math>. | ||
+ | |||
+ | Таким образом, в параллельном варианте, в отличие от последовательного, вычисления квадратных корней и делений будут определять довольно значительную долю требуемого времени. При реализации на конкретных архитектурах наличие в отдельных ярусах [[глоссарий#Ярусно-параллельная форма графа алгоритма|ЯПФ]] отдельных вычислений квадратных корней может породить и другие проблемы. Например, при реализации на ПЛИСах остальные вычисления (деления и тем более умножения и сложения/вычитания) могут быть конвейеризованы, что даёт экономию и по ресурсам на программируемых платах; вычисления же квадратных корней из-за их изолированности приведут к занятию ресурсов на платах, которые будут простаивать большую часть времени. Таким образом, общая экономия в 2 раза, из-за которой метод Холецкого предпочитают в случае симметричных задач тому же методу Гаусса, в параллельном случае уже имеет место вовсе не по всем ресурсам, и главное - не по требуемому времени. | ||
+ | |||
+ | При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности, а в параллельном варианте это означает, что практически все промежуточные вычисления для выполнения метода Холецкого в режиме накопления должны быть двойной точности. В отличие от последовательного варианта это означает увеличение требуемой памяти почти в 2 раза. | ||
+ | |||
+ | При классификации по высоте ЯПФ, таким образом, метод Холецкого относится к алгоритмам со сложностью <math>O(n)</math>. При классификации по ширине ЯПФ его сложность будет <math>O(n^2)</math>. | ||
+ | |||
+ | === Входные и выходные данные алгоритма === | ||
+ | |||
+ | '''Входные данные''': плотная матрица <math>A</math> (элементы <math>a_{ij}</math>). | ||
+ | Дополнительные ограничения: | ||
+ | * <math>A</math> – симметрическая матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>. | ||
+ | * <math>A</math> – положительно определённая матрица, т. е. для любых ненулевых векторов <math>\vec{x}</math> выполняется <math>\vec{x}^T A \vec{x} > 0</math>. | ||
+ | |||
+ | '''Объём входных данных''': <math>\frac{n (n + 1)}{2}</math> (в силу симметричности достаточно хранить только диагональ и над/поддиагональные элементы). В разных реализациях эта экономия хранения может быть выполнена разным образом. Например, в библиотеке, реализованной в НИВЦ МГУ, матрица A хранилась в одномерном массиве длины <math>\frac{n (n + 1)}{2}</math> по строкам своего нижнего треугольника. | ||
+ | |||
+ | '''Выходные данные''': нижняя треугольная матрица <math>L</math> (элементы <math>l_{ij}</math>). | ||
+ | |||
+ | '''Объём выходных данных''': <math>\frac{n (n + 1)}{2}</math> (в силу треугольности достаточно хранить только ненулевые элементы). В разных реализациях эта экономия хранения может быть выполнена разным образом. Например, в той же библиотеке, созданной в НИВЦ МГУ, матрица <math>L</math> хранилась в одномерном массиве длины <math>\frac{n (n + 1)}{2}</math> по строкам своей нижней части. | ||
+ | |||
+ | === Свойства алгоритма === | ||
+ | |||
+ | Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является ''квадратичным'' (отношение кубической к линейной). | ||
+ | |||
+ | При этом вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных – всего лишь ''линейна''. | ||
+ | |||
+ | При этом алгоритм почти полностью детерминирован, это гарантируется теоремой о единственности разложения. Использование другого порядка выполнения ассоциативных операций может привести к накоплению ошибок округления, однако это влияние в используемых вариантах алгоритма не так велико, как, скажем, отказ от использования режима накопления. | ||
+ | |||
+ | Дуги информационного графа, исходящие из вершин, соответствующих операциям квадратного корня и деления, образуют пучки т. н. рассылок линейной мощности (то есть степень исхода этих вершин и мощность работы с этими данными — линейная функция от порядка матрицы и координат этих вершин). При этом естественно наличие в этих пучках «длинных» дуг. Остальные дуги локальны. | ||
+ | |||
+ | Наиболее известной является компактная укладка графа — его проекция на треугольник матрицы, который перевычисляется укладываемыми операциями. При этом «длинные» дуги можно убрать, заменив более дальнюю пересылку комбинацией нескольких ближних (к соседям). | ||
+ | |||
+ | [[Глоссарий#Эквивалентное возмущение|Эквивалентное возмущение]] <math>M</math> у метода Холецкого всего вдвое больше, чем возмущение <math>\delta A</math>, вносимое в матрицу при вводе чисел в компьютер: | ||
+ | <math> | ||
+ | ||M||_{E} \leq 2||\delta A||_{E} | ||
+ | </math> | ||
+ | |||
+ | Это явление обусловлено положительной определённостью матрицы. Среди всех используемых разложений матриц это наименьшее из эквивалентных возмущений. | ||
+ | |||
+ | == Литература == | ||
+ | |||
+ | <references \> |
Версия 11:04, 17 сентября 2015
Содержание
- 1 Разложение Холецкого (метод квадратного корня)
- 1.1 Свойства и структура алгоритма
- 1.1.1 Общее описание алгоритма
- 1.1.2 Математическое описание алгоритма
- 1.1.3 Вычислительное ядро алгоритма
- 1.1.4 Макроструктура алгоритма
- 1.1.5 Схема реализации последовательного алгоритма
- 1.1.6 Последовательная сложность алгоритма
- 1.1.7 Информационный граф
- 1.1.8 Ресурс параллелизма алгоритма
- 1.1.9 Входные и выходные данные алгоритма
- 1.1.10 Свойства алгоритма
- 1.2 Литература
- 1.1 Свойства и структура алгоритма
1 Разложение Холецкого (метод квадратного корня)
1.1 Свойства и структура алгоритма
1.1.1 Общее описание алгоритма
Разложение Холецкого впервые предложено французским офицером и математиком Андре-Луи Холецким в конце Первой Мировой войны, незадолго до его гибели в бою в августе 1918 г. Идея этого разложения была опубликована в 1924 г. его сослуживцем[1]. Потом оно было использовано поляком Т. Банашевичем[2][3] в 1938 г. В советской математической литературе называется также методом квадратного корня[4][5][6]; название связано с характерными операциями, отсутствующими в родственном разложении Гаусса.
Первоначально разложение Холецкого использовалось исключительно для плотных симметричных положительно определенных матриц. В настоящее время его использование гораздо шире. Оно может быть применено также, например, к эрмитовым матрицам. Для повышения производительности вычислений часто применяется блочная версия разложения.
Для разреженных матриц разложение Холецкого также широко применяется в качестве основного этапа прямого метода решения линейных систем. В этом случае используют специальные упорядочивания для уменьшения ширины профиля исключения, а следовательно и уменьшения количества арифметических операций. Другие упорядочивания используются для выделения независимых блоков вычислений при работе на системах с параллельной организацией.
Варианты разложения Холецкого нашли успешные применения и в итерационных методах для построения переобусловливателей разреженных симметричных положительно определенных матриц. В неполном треугольном разложении («по позициям») элементы переобусловливателя вычисляются только в заранее заданных позициях, например, в позициях ненулевых элементов исходной матрицы (так называемое разложение IC0). Для получения же более точного разложения применяется приближение, в котором фильтрация малых элементов производится «по значениям». В зависимости от используемого порога фильтрации можно получить более точное, но и более заполненное разложение. Существует и алгоритм разложения второго порядка точности[7]. В нём при таком же заполнении множителей разложения удается улучшить точность. Для такого разложения в параллельном режиме используется специальный вариант аддитивного переобуславливания на основе разложения второго порядка[8].
На этой странице представлено исходное разложение Холецкого с новых позиций нашего суперкомпьютерного века. Приведено описание конкретной версии разложения Холецкого — для плотных вещественных симметричных положительно определённых матриц, но структура для ряда других версий, например, для комплексного случая, почти такая же, различия состоят в замене большинства вещественных операций на комплексные. Список остальных основных вариантов разложения можно посмотреть на странице Метод Холецкого (нахождение симметричного треугольного разложения).
Используется для разложения положительно определённых эрмитовых (в вещественном случае - симметрических) матриц в виде [math]A = L L^*[/math] ([math]L[/math] — нижняя треугольная матрица) или в виде [math]A = U^* U[/math] ([math]U[/math] — верхняя треугольная матрица ; эти разложения совершенно эквивалентны друг другу по вычислениям и различаются только способом представления данных). Он заключается в реализации формул для элементов матрицы [math]L[/math], получающихся из вышеприведённого равенства единственным образом. Получило широкое распространение благодаря следующим особенностям.
1.1.1.1 Симметричность и положительная определённость матрицы
Симметричность матрицы позволяет хранить и вычислять только чуть больше половины её элементов, что почти вдвое экономит как необходимые для вычислений объёмы памяти, так и количество операций в сравнении с, например, разложением по методу Гаусса. При этом альтернативное (без вычисления квадратных корней) LU-разложение, использующее симметрию матрицы, всё же несколько быстрее метода Холецкого (не использует извлечение квадратных корней), но требует хранения всей матрицы. Благодаря тому, что разлагаемая матрица не только симметрична, но и положительно определена, её LU-разложения, в том числе и разложение методом Холецкого, имеют наименьшее эквивалентное возмущение из всех известных разложений матриц.
1.1.1.2 Режим накопления
Алгоритм позволяет использовать так называемый режим накопления, обусловленный тем, что значительную часть вычислений составляют вычисления скалярных произведений.
1.1.2 Математическое описание алгоритма
Исходные данные: положительно определённая симметрическая матрица [math]A[/math] (элементы [math]a_{ij}[/math]).
Вычисляемые данные: нижняя треугольная матрица [math]L[/math] (элементы [math]l_{ij}[/math]).
Формулы метода:
- [math] \begin{align} l_{11} & = \sqrt{a_{11}}, \\ l_{j1} & = \frac{a_{j1}}{l_{11}}, \quad j \in [2, n], \\ l_{ii} & = \sqrt{a_{ii} - \sum_{p = 1}^{i - 1} l_{ip}^2}, \quad i \in [2, n], \\ l_{ji} & = \left (a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp} \right ) / l_{ii}, \quad i \in [2, n - 1], j \in [i + 1, n]. \end{align} [/math]
Существует также блочная версия метода, однако в данном описании разобран только точечный метод.
В ряде реализаций деление на диагональный элемент выполняется в два этапа: вычисление [math]\frac{1}{l_{ii}}[/math] и затем умножение на него всех (видоизменённых) [math]a_{ji}[/math] . Здесь мы этот вариант алгоритма не рассматриваем. Заметим только, что он имеет худшие параллельные характеристики, чем представленный.
1.1.3 Вычислительное ядро алгоритма
Вычислительное ядро последовательной версии метода Холецкого можно составить из множественных (всего их [math]\frac{n (n - 1)}{2}[/math]) вычислений скалярных произведений строк матрицы:
- [math]\sum_{p = 1}^{i - 1} l_{ip} l_{jp}[/math]
в режиме накопления или без него, в зависимости от требований задачи. Во многих последовательных реализациях упомянутый способ представления не используется. Дело в том, что в них вычисление сумм типа
- [math]a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp}[/math]
в которых и встречаются скалярные произведения, ведутся не в порядке «вычислили скалярное произведение, а потом вычли его из элемента», а путём вычитания из элемента покомпонентных произведений, являющихся частями скалярных произведений. Поэтому следует считать вычислительным ядром метода не вычисления скалярных произведений, а вычисления выражений
- [math]a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp}[/math]
в режиме накопления или без него.
Тем не менее, в популярных зарубежных реализациях точечного метода Холецкого, в частности, в библиотеках LINPACK и LAPACK, основанных на BLAS, используются именно вычисления скалярных произведений в виде вызова соответствующих подпрограмм BLAS (конкретно — функции SDOT). На последовательном уровне это влечёт за собой дополнительную операцию суммирования на каждый из [math]\frac{n (n + 1)}{2}[/math] вычисляемый элемент матрицы [math]L[/math] и некоторое замедление работы программы (о других следствиях рассказано ниже в разделе «Существующие реализации алгоритма»). Поэтому в данных вариантах ядром метода Холецкого будут вычисления этих скалярных произведений.
1.1.4 Макроструктура алгоритма
Как записано и в описании ядра алгоритма, основную часть метода составляют множественные (всего [math]\frac{n (n - 1)}{2}[/math]) вычисления сумм
- [math]a_{ji}-\sum_{p=1}^{i-1}l_{ip} l_{jp}[/math]
в режиме накопления или без него.
1.1.5 Схема реализации последовательного алгоритма
Последовательность исполнения метода следующая:
1. [math]l_{11}= \sqrt{a_{11}}[/math]
2. [math]l_{j1}= \frac{a_{j1}}{l_{11}}[/math] (при [math]j[/math] от [math]2[/math] до [math]n[/math]).
Далее для всех [math]i[/math] от [math]2[/math] до [math]n[/math] по нарастанию выполняются
3. [math]l_{ii} = \sqrt{a_{ii} - \sum_{p = 1}^{i - 1} l_{ip}^2}[/math] и
4. (кроме [math]i = n[/math]): [math]l_{ji} = \left (a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp} \right ) / l_{ii}[/math] (для всех [math]j[/math] от [math]i + 1[/math] до [math]n[/math]).
После этого (если [math]i \lt n[/math]) происходит переход к шагу 3 с бо́льшим [math]i[/math].
Особо отметим, что вычисления сумм вида [math]a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp}[/math] в обеих формулах производят в режиме накопления вычитанием из [math]a_{ji}[/math] произведений [math]l_{ip} l_{jp}[/math] для [math]p[/math] от [math]1[/math] до [math]i - 1[/math], c нарастанием [math]p[/math].
1.1.6 Последовательная сложность алгоритма
Для разложения матрицы порядка n методом Холецкого в последовательном (наиболее быстром) варианте требуется:
- [math]n[/math] вычислений квадратного корня,
- [math]\frac{n(n-1)}{2}[/math] делений,
- [math]\frac{n^3-n}{6}[/math] сложений (вычитаний),
- [math]\frac{n^3-n}{6}[/math] умножений.
Умножения и сложения (вычитания) составляют основную часть алгоритма.
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности (или использования функции вроде DPROD в Фортране), что ещё больше увеличивает долю умножений и сложений/вычитаний во времени, требуемом для выполнения метода Холецкого.
При классификации по последовательной сложности, таким образом, метод Холецкого относится к алгоритмам с кубической сложностью.
1.1.7 Информационный граф
Опишем граф алгоритма[9][10][11] как аналитически, так и в виде рисунка.
Граф алгоритма состоит из трёх групп вершин, расположенных в целочисленных узлах трёх областей разной размерности.
Первая группа вершин расположена в одномерной области, соответствующая ей операция вычисляет функцию SQRT. Единственная координата каждой из вершин [math]i[/math] меняется в диапазоне от [math]1[/math] до [math]n[/math], принимая все целочисленные значения.
Аргумент этой функции
- при [math]i = 1[/math] — элемент входных данных, а именно [math]a_{11}[/math];
- при [math]i \gt 1[/math] — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами [math]i - 1[/math], [math]i[/math], [math]i - 1[/math].
Результат срабатывания операции является выходным данным [math]l_{ii}[/math].
Вторая группа вершин расположена в двумерной области, соответствующая ей операция [math]a / b[/math]. Естественно введённые координаты области таковы:
- [math]i[/math] — меняется в диапазоне от [math]1[/math] до [math]n-1[/math], принимая все целочисленные значения;
- [math]j[/math] — меняется в диапазоне от [math]i+1[/math] до [math]n[/math], принимая все целочисленные значения.
Аргументы операции следующие:
- [math]a[/math]:
- при [math]i = 1[/math] — элементы входных данных, а именно [math]a_{j1}[/math];
- при [math]i \gt 1[/math] — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами [math]i - 1, j, i - 1[/math];
- [math]b[/math] — результат срабатывания операции, соответствующей вершине из первой группы, с координатой [math]i[/math].
Результат срабатывания операции является выходным данным [math]l_{ji}[/math].
Третья группа вершин расположена в трёхмерной области, соответствующая ей операция [math]a - b * c[/math]. Естественно введённые координаты области таковы:
- [math]i[/math] — меняется в диапазоне от [math]2[/math] до [math]n[/math], принимая все целочисленные значения;
- [math]j[/math] — меняется в диапазоне от [math]i[/math] до [math]n[/math], принимая все целочисленные значения;
- [math]p[/math] — меняется в диапазоне от [math]1[/math] до [math]i - 1[/math], принимая все целочисленные значения.
Аргументы операции следующие:
- [math]a[/math]:
- при [math]p = 1[/math] элемент входных данных [math]a_{ji}[/math];
- при [math]p \gt 1[/math] — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами [math]i, j, p - 1[/math];
- [math]b[/math] — результат срабатывания операции, соответствующей вершине из второй группы, с координатами [math]p, i[/math];
- [math]c[/math] — результат срабатывания операции, соответствующей вершине из второй группы, с координатами [math]p, j[/math];
Результат срабатывания операции является промежуточным данным алгоритма.
Описанный граф можно посмотреть на рис.1 и рис.2, выполненных для случая [math]n = 4[/math]. Здесь вершины первой группы обозначены жёлтым цветом и буквосочетанием sq, вершины второй — зелёным цветом и знаком деления, третьей — красным цветом и буквой f. Вершины, соответствующие операциям, производящим выходные данные алгоритма, выполнены более крупно. Дублирующие друг друга дуги даны как одна. На рис.1 показан граф алгоритма согласно классическому определению , на рис.2 к графу алгоритма добавлены вершины , соответствующие входным (обозначены синим цветом) и выходным (обозначены розовым цветом) данным.
1.1.8 Ресурс параллелизма алгоритма
Для разложения матрицы порядка [math]n[/math] методом Холецкого в параллельном варианте требуется последовательно выполнить следующие ярусы:
- [math]n[/math] ярусов с вычислением квадратного корня (единичные вычисления в каждом из ярусов),
- [math]n - 1[/math] ярус делений (в каждом из ярусов линейное количество делений, в зависимости от яруса — от [math]1[/math] до [math]n - 1[/math]),
- по [math]n - 1[/math] ярусов умножений и сложений/вычитаний (в каждом из ярусов — квадратичное количество операций, от [math]1[/math] до [math]\frac{n^2 - n}{2}[/math].
Таким образом, в параллельном варианте, в отличие от последовательного, вычисления квадратных корней и делений будут определять довольно значительную долю требуемого времени. При реализации на конкретных архитектурах наличие в отдельных ярусах ЯПФ отдельных вычислений квадратных корней может породить и другие проблемы. Например, при реализации на ПЛИСах остальные вычисления (деления и тем более умножения и сложения/вычитания) могут быть конвейеризованы, что даёт экономию и по ресурсам на программируемых платах; вычисления же квадратных корней из-за их изолированности приведут к занятию ресурсов на платах, которые будут простаивать большую часть времени. Таким образом, общая экономия в 2 раза, из-за которой метод Холецкого предпочитают в случае симметричных задач тому же методу Гаусса, в параллельном случае уже имеет место вовсе не по всем ресурсам, и главное - не по требуемому времени.
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности, а в параллельном варианте это означает, что практически все промежуточные вычисления для выполнения метода Холецкого в режиме накопления должны быть двойной точности. В отличие от последовательного варианта это означает увеличение требуемой памяти почти в 2 раза.
При классификации по высоте ЯПФ, таким образом, метод Холецкого относится к алгоритмам со сложностью [math]O(n)[/math]. При классификации по ширине ЯПФ его сложность будет [math]O(n^2)[/math].
1.1.9 Входные и выходные данные алгоритма
Входные данные: плотная матрица [math]A[/math] (элементы [math]a_{ij}[/math]). Дополнительные ограничения:
- [math]A[/math] – симметрическая матрица, т. е. [math]a_{ij}= a_{ji}, i, j = 1, \ldots, n[/math].
- [math]A[/math] – положительно определённая матрица, т. е. для любых ненулевых векторов [math]\vec{x}[/math] выполняется [math]\vec{x}^T A \vec{x} \gt 0[/math].
Объём входных данных: [math]\frac{n (n + 1)}{2}[/math] (в силу симметричности достаточно хранить только диагональ и над/поддиагональные элементы). В разных реализациях эта экономия хранения может быть выполнена разным образом. Например, в библиотеке, реализованной в НИВЦ МГУ, матрица A хранилась в одномерном массиве длины [math]\frac{n (n + 1)}{2}[/math] по строкам своего нижнего треугольника.
Выходные данные: нижняя треугольная матрица [math]L[/math] (элементы [math]l_{ij}[/math]).
Объём выходных данных: [math]\frac{n (n + 1)}{2}[/math] (в силу треугольности достаточно хранить только ненулевые элементы). В разных реализациях эта экономия хранения может быть выполнена разным образом. Например, в той же библиотеке, созданной в НИВЦ МГУ, матрица [math]L[/math] хранилась в одномерном массиве длины [math]\frac{n (n + 1)}{2}[/math] по строкам своей нижней части.
1.1.10 Свойства алгоритма
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является квадратичным (отношение кубической к линейной).
При этом вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных – всего лишь линейна.
При этом алгоритм почти полностью детерминирован, это гарантируется теоремой о единственности разложения. Использование другого порядка выполнения ассоциативных операций может привести к накоплению ошибок округления, однако это влияние в используемых вариантах алгоритма не так велико, как, скажем, отказ от использования режима накопления.
Дуги информационного графа, исходящие из вершин, соответствующих операциям квадратного корня и деления, образуют пучки т. н. рассылок линейной мощности (то есть степень исхода этих вершин и мощность работы с этими данными — линейная функция от порядка матрицы и координат этих вершин). При этом естественно наличие в этих пучках «длинных» дуг. Остальные дуги локальны.
Наиболее известной является компактная укладка графа — его проекция на треугольник матрицы, который перевычисляется укладываемыми операциями. При этом «длинные» дуги можно убрать, заменив более дальнюю пересылку комбинацией нескольких ближних (к соседям).
Эквивалентное возмущение [math]M[/math] у метода Холецкого всего вдвое больше, чем возмущение [math]\delta A[/math], вносимое в матрицу при вводе чисел в компьютер: [math] ||M||_{E} \leq 2||\delta A||_{E} [/math]
Это явление обусловлено положительной определённостью матрицы. Среди всех используемых разложений матриц это наименьшее из эквивалентных возмущений.
1.2 Литература
<references \>
- ↑ Commandant Benoit, Note sur une méthode de résolution des équations normales provenant de l'application de la méthode des moindres carrés à un système d'équations linéaires en nombre inférieur à celui des inconnues (Procédé du Commandant Cholesky), Bulletin Géodésique 2 (1924), 67-77.
- ↑ Banachiewicz T. Principles d'une nouvelle technique de la méthode des moindres carrês. Bull. Intern. Acad. Polon. Sci. A., 1938, 134-135.
- ↑ Banachiewicz T. Méthode de résoltution numérique des équations linéaires, du calcul des déterminants et des inverses et de réduction des formes quardatiques. Bull. Intern. Acad. Polon. Sci. A., 1938, 393-401.
- ↑ Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977.
- ↑ Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. М.: Наука, 1984.
- ↑ Фаддеев Д.К., Фаддева В.Н. Вычислительные основы линейной алгебры. М.-Л.: Физматгиз, 1963.
- ↑ Kaporin I.E. High quality preconditioning of a general symmetric positive definite matrix based on its UTU + UTR + RTU-decomposition. Numer. Lin. Algebra Appl. (1998) Vol. 5, No. 6, 483-509.
- ↑ Капорин И.Е., Коньшин И.Н. Параллельное решение симметричных положительно-определенных систем на основе перекрывающегося разбиения на блоки. Ж. вычисл. матем. и матем. физ., 2001, Т, 41, N. 4, C. 515–528.
- ↑ Воеводин В.В. Математические основы параллельных вычислений// М.: Изд. Моск. ун-та, 1991. 345 с.
- ↑ Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. – СПб.: БХВ - Петербург, 2002. – 608 с.
- ↑ Фролов А.В.. Принципы построения и описание языка Сигма. Препринт ОВМ АН N 236. М.: ОВМ АН СССР, 1989.