|
|
| Строка 1: |
Строка 1: |
| − | {{algorithm
| + | [[Приложение 1]] |
| − | | name = Разложение Холецкого
| |
| − | | serial_complexity = <math>O(n^3)</math>
| |
| − | | pf_height = <math>O(n)</math>
| |
| − | | pf_width = <math>O(n^2)</math>
| |
| − | | input_data = <math>\frac{n (n + 1)}{2}</math>
| |
| − | | output_data = <math>\frac{n (n + 1)}{2}</math>
| |
| − | }}
| |
| − | | |
| − | == Свойства и структура алгоритма ==
| |
| − | | |
| − | === Общее описание алгоритма ===
| |
| − | | |
| − | '''Разложение Холецкого''' впервые предложено французским офицером и математиком Андре-Луи Холецким в конце Первой Мировой войны, незадолго до его гибели в бою в августе 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 \>
| |
