https://algowiki-project.org/w/ru/api.php?action=feedcontributions&user=Oleg+Arushanyan&feedformat=atomАлговики - Вклад участника [ru]2024-03-28T17:32:46ZВклад участникаMediaWiki 1.34.0https://algowiki-project.org/w/ru/index.php?title=%D0%9E%D0%B1%D1%80%D0%B0%D1%82%D0%BD%D0%B0%D1%8F_%D0%BF%D0%BE%D0%B4%D1%81%D1%82%D0%B0%D0%BD%D0%BE%D0%B2%D0%BA%D0%B0_(%D0%B2%D0%B5%D1%89%D0%B5%D1%81%D1%82%D0%B2%D0%B5%D0%BD%D0%BD%D1%8B%D0%B9_%D0%B2%D0%B0%D1%80%D0%B8%D0%B0%D0%BD%D1%82)&diff=2002Обратная подстановка (вещественный вариант)2015-06-28T14:54:13Z<p>Oleg Arushanyan: /* Общее описание алгоритма */</p>
<hr />
<div>Основные авторы описания: [[Участник:Frolov|А.В.Фролов]], [[Участник:VadimVV|Вад.В.Воеводин]] ([[#Описание локальности данных и вычислений|раздел 2.2]])<br />
<br />
== Описание свойств и структуры алгоритма ==<br />
<br />
=== Общее описание алгоритма ===<br />
<br />
'''Обратная подстановка''' - решение ''системы линейных алгебраических уравнений'' ('''СЛАУ''') <math>Ux = y</math> с верхней треугольной матрицей <math>U</math>. Матрица <math>U</math> может быть одной из составляющих матрицы <math>A</math> в каких-либо разложениях и получается либо из <math>LU</math>-разложения последней каким-либо из многочисленных способов (например, простое разложение Гаусса, разложение Гаусса с выбором ведущего элемента, компактная схема Гаусса, [[Метод Холецкого (квадратного корня), точечный вещественный вариант|разложение Холецкого]] и др.), либо из других (например из QR-разложения). В силу треугольности <math>U</math> решение СЛАУ является одной из модификаций общего метода подстановки и записывается простыми формулами.<br />
<br />
В [1] методом '''обратной подстановки''' назван также и метод решения СЛАУ с ''нижней треугольной матрицей''. Там же отмечено, что в литературе иногда под ''обратной подстановкой'' имеют в виду, как и здесь, только решения СЛАУ с ''верхней треугольной матрицей'', а решение ''нижних'' треугольных систем называют [[Прямая подстановка (вещественный вариант)|прямой подстановкой]]. Такой же системы названий будем придерживаться и здесь, во избежание одноимённого названия разных алгоритмов. Кроме того, '''обратная подстановка''', представленная здесь, одновременно может быть частью '''метода Гаусса для решения СЛАУ''', а именно - его '''обратным ходом''', чего нельзя сказать про [[Прямая подстановка (вещественный вариант)|прямую подстановку]].<br />
<br />
Существует метод со сходным названием - [[Обратная подстановка с нормировкой]]. При том, что он решает, по существу, ту же задачу, что и простая '''обратная подстановка''', его схема несколько сложнее. Это связано со специальными мерами по уменьшению влияния ошибок округления на результат. [[Обратная подстановка с нормировкой]] на данной странице не рассматривается.<br />
<br />
=== Математическое описание ===<br />
<br />
Исходные данные: верхняя треугольная матрица <math>U</math> (элементы <math>u_{ij}</math>), вектор правой части <math>y</math> (элементы <math>y_{i}</math>).<br />
<br />
Вычисляемые данные: вектор решения <math>x</math> (элементы <math>x_{i}</math>).<br />
<br />
Формулы метода:<br />
:<math><br />
\begin{align}<br />
x_{n} & = y_{n}/u_{nn} \\<br />
x_{i} & = \left (y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} \right ) / u_{ii}, \quad i \in [1, n - 1].<br />
\end{align}<br />
</math><br />
<br />
Существует также блочная версия метода, однако в данном описании разобран только точечный метод.<br />
<br />
=== Вычислительное ядро алгоритма ===<br />
<br />
Вычислительное ядро обратной подстановки можно составить из множественных (всего их <math>n-1</math>) вычислений скалярных произведений подстрок матрицы <math>U</math> на уже вычисленную часть вектора <math>x</math>:<br />
<br />
:<math> \sum_{j = i+1}^{n} u_{ij} x_{j} </math> <br />
<br />
в режиме накопления или без него, в зависимости от требований задачи, с их последующим вычитанием из компоненты вектора <math>y</math> и деления на диагональный элемент матрицы <math>U</math>. В отечественных реализациях, даже в последовательных, упомянутый способ представления не используется. Дело в том, что даже в этих реализациях метода вычисление сумм типа <br />
<br />
:<math> y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} </math><br />
<br />
в которых и встречаются скалярные произведения, ведутся не в порядке «вычислили скалярное произведение, а потом вычли его из элемента», а путём вычитания из элемента покомпонентных произведений, являющихся частями скалярных произведений. Поэтому следует считать вычислительным ядром метода не вычисления скалярных произведений, а вычисления выражений<br />
<br />
:<math> y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} </math><br />
<br />
в режиме накопления или без него, в зависимости от требований задачи, плюс деления результатов этих вычислений на диагональные элементы матрицы.<br />
<br />
=== Макроструктура алгоритма ===<br />
<br />
Как уже записано в [[#Вычислительное ядро алгоритма|описании ядра алгоритма]], основную часть метода обратной подстановки составляют множественные (всего <math>n-1</math>) вычисления сумм<br />
<br />
:<math>y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} </math><br />
<br />
в режиме накопления или без него, плюс деления результатов этих вычислений на диагональные элементы матрицы.<br />
<br />
=== Описание схемы реализации последовательного алгоритма ===<br />
<br />
Последовательность исполнения такова:<br />
<br />
1. <math>x_{n} = y_{n}/u_{nn}</math><br />
<br />
Далее для всех <math>i</math> от <math>n-1</math> до <math>1</math> ''по убыванию'' выполняются<br />
<br />
2. <math>x_{i} = \left (y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} \right ) / u_{ii}</math><br />
<br />
Особо отметим, что вычисления сумм вида <math>y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j}</math> производят в режиме накопления вычитанием из <math>y_{i}</math> произведений <math>u_{ij} x_{j}</math> для <math>j</math> от <math>n</math> до <math>i + 1</math>, '''''c убыванием''''' <math>j</math>. '''''Другие порядки выполнения суммирования приводят к резкому ухудшению параллельных свойств алгоритма''''', хотя, к сожалению, остаются кое-где в литературе и пакетах программ. В качестве примера такого порядка можно привести фрагмент программы из [2], где обратная подстановка является обратным ходом в методе Гаусса, а возрастание индекса суммирования связано, в основном, с ограничениями используемого авторами книги старого диалекта Фортрана.<br />
<br />
=== Последовательная сложность алгоритма ===<br />
<br />
Для обратной подстановки порядка n в последовательном (наиболее быстром) варианте требуется:<br />
<br />
* <math>n</math> делений,<br />
* <math>\frac{n^2-n}{2}</math> сложений (вычитаний),<br />
* <math>\frac{n^2-n}{2}</math> умножений. <br />
Умножения и сложения (вычитания) — ''основная часть алгоритма''.<br />
<br />
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности (или использования функции вроде DPROD в Фортране).<br />
<br />
При классификации по последовательной сложности, таким образом, метод обратной подстановки относится к алгоритмам ''со сложностью'' <math>O(n^2)</math>.<br />
<br />
=== Информационный граф ===<br />
<br />
Опишем [[глоссарий#Граф алгоритма|граф алгоритма]] как аналитически, так и в виде рисунка. <br />
<br />
Граф алгоритма обратной подстановки состоит из двух групп вершин, расположенных в целочисленных узлах двух областей разной размерности.<br />
<br />
[[Файл:DirectU.png|450px|thumb|left|Обратная подстановка]]<br />
<br />
'''Первая''' группа вершин расположена в одномерной области, соответствующая ей операция вычисляет функцию деления. <br />
Естественно введённая единственная координата каждой из вершин <math>i</math> меняется в диапазоне от <math>n</math> до <math>1</math>, принимая все целочисленные значения.<br />
<br />
Делимое в этой операции: <br />
<br />
* при <math>i = n</math> — элемент ''входных данных'', а именно <math>y_{n}</math>; <br />
* при <math>i < n</math> — результат срабатывания операции, соответствующей вершине из второй группы, с координатами <math>i</math>, <math>i+1</math>.<br />
<br />
Делитель для этой операции - элемент ''входных данных'', а именно <math>u_{nn}</math>.<br />
<br />
Результат срабатывания операции является ''выходным данным'' <math>x_{i}</math>.<br />
<br />
'''Вторая''' группа вершин расположена в двумерной области, соответствующая ей операция <math>a-bc</math>. <br />
Естественно введённые координаты области таковы: <br />
* <math>i</math> — меняется в диапазоне от <math>n-1</math> до <math>1</math>, принимая все целочисленные значения;<br />
* <math>j</math> — меняется в диапазоне от <math>n</math> до <math>i+1</math>, принимая все целочисленные значения.<br />
<br />
Аргументы операции следующие:<br />
*<math>a</math>:<br />
** при <math>j = n</math> элемент входных данных <math>y_{i}</math>; <br />
** при <math>j < n</math> — результат срабатывания операции, соответствующей вершине из второй группы, с координатами <math>i, j+1</math>;<br />
*<math>b</math> — элемент ''входных данных'', а именно <math>u_{ij}</math>;<br />
*<math>c</math> — результат срабатывания операции, соответствующей вершине из первой группы, с координатой <math>j</math>;<br />
<br />
Результат срабатывания операции является ''промежуточным данным'' алгоритма.<br />
<br />
Описанный граф можно посмотреть на рисунке, выполненном для случая <math>n = 5</math>. Здесь вершины первой группы обозначены ;`жёлтым цветом и знаком деления, вершины второй — зелёным цветом и буквой f. Изображена подача только входных данных из вектора <math>y</math>, подача элементов матрицы <math>U</math>, идущая во все вершины, на рисунке не представлена.<br />
<br />
=== Описание ресурса параллелизма алгоритма ===<br />
<br />
Для обратной подстановки порядка n в параллельном варианте требуется последовательно выполнить следующие ярусы:<br />
* <math>n</math> ярусов делений (в каждом из ярусов одно деление),<br />
* по <math>n - 1</math> ярусов умножений и сложений/вычитаний (в каждом из ярусов — линейное количество операций, от <math>1</math> до <math>n-1</math>.<br />
<br />
Таким образом, в параллельном варианте, в отличие от последовательного, вычисления делений будут определять довольно значительную долю требуемого времени. При реализации на конкретных архитектурах наличие в отдельных ярусах [[глоссарий#Ярусно-параллельная форма графа алгоритма|ЯПФ]] отдельных делений может породить и другие проблемы. Например, при реализации метода обратной подстановки на ПЛИСах остальные вычисления (умножения и сложения/вычитания) могут быть конвейеризованы, что даёт экономию и по ресурсам на программируемых платах; деления из-за их изолированности приведут к занятию ресурсов на платах, которые будут простаивать большую часть времени. <br />
<br />
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности, а в параллельном варианте это означает, что практически все промежуточные вычисления для выполнения алгоритма в режиме накопления должны быть двойной точности. В отличие от последовательного варианта это означает некоторое увеличение требуемой памяти.<br />
<br />
При классификации по высоте ЯПФ, таким образом, метод обратной подстановки относится к алгоритмам ''с линейной сложностью''. При классификации по ширине ЯПФ его сложность также будет ''линейной''.<br />
<br />
=== Описание входных и выходных данных ===<br />
<br />
'''Входные данные''': верхняя треугольная матрица <math>U</math> (элементы <math>u_{ij}</math>), вектор правой части <math>y</math> (элементы <math>y_{i}</math>).<br />
<br />
'''Объём входных данных''': :<math>\frac{n (n + 3)}{2}</math> (в силу треугольности достаточно хранить только ненулевые элементы матрицы <math>U</math>). <br />
<br />
'''Выходные данные''': вектор решения <math>x</math> (элементы <math>x_{i}</math>).<br />
<br />
'''Объём выходных данных''': :<math>n~.</math><br />
<br />
=== Свойства алгоритма ===<br />
<br />
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является ''линейным'' (отношение квадратической к линейной). <br />
<br />
При этом вычислительная мощность алгоритма обратной подстановки, как отношение числа операций к суммарному объему входных и выходных данных – всего лишь ''константа''.<br />
<br />
При этом алгоритм обратной подстановки полностью детерминирован. Использование другого порядка выполнения ассоциативных операций в данной версии нами не рассматривается, поскольку в корне меняет структуру алгоритма и меняет сложность с параллельной на квадратичную.<br />
<br />
Наличие линейного количества ярусов ЯПФ, состоящих из одного-единственного деления, потенциально замедляющее параллельные реализации алгоритма, является его характерным "узким местом", особенно в сравнении со схожей по решаемой математической задаче [[Прямая подстановка (вещественный вариант)|прямой подстановке]], где диагональные элементы единичны. В связи с этим для решения СЛАУ предпочтительны такие разложения, содержащие треугольные матрицы, где в треугольных матрицах диагональные элементы единичны. В тех же случаях, когда получаются неособенные треугольные матрицы, их желательно предварительно, до решения СЛАУ с ними, преобразовать в произведение диагональной и треугольной с единичными диагональными элементами.<br />
<br />
У алгоритма обратной подстановки существует несколько блочных вариантов. Граф некоторых из них совпадает с графом точечного варианта, различия связаны в основном с порядком прохождения основных циклов алгоритма, а именно - с их развёртыванием и перестановкой. Эти приёмы могут помочь в оптимизации обменов на конкретных вычислительных системах.<br />
<br />
== Программная реализация ==<br />
<br />
=== Особенности реализации последовательного алгоритма ===<br />
<br />
В простейшем варианте метод обратной подстановки на Фортране можно записать так:<br />
<br />
<source lang="fortran"><br />
X(N) = Y(N) / U (N, N)<br />
DO I = N-1, 1, -1<br />
S = Y(I)<br />
DO J = N, I+1, -1<br />
S = S - DPROD(U(I,J), X(J))<br />
END DO<br />
X(I) = S / U(I,I) <br />
END DO<br />
</source><br />
При этом для реализации режима накопления переменная <math>S</math> должна быть двойной точности.<br />
<br />
=== Описание локальности данных и вычислений ===<br />
==== Описание локальности реализации алгоритма ====<br />
===== Описание структуры обращений в память и качественная оценка локальности =====<br />
<br />
[[file:gauss back 1.PNG|thumb|center|700px|Рисунок 12.1. Обратный ход метода Гаусса решения СЛАУ. Общий профиль обращений в память]]<br />
<br />
На рисунке 12.1 представлен профиль обращений в память для обратного хода метода Гаусса решения СЛАУ. Хорошо видно, что профиль состоит из двух этапов, идущих друг за другом (граница между ними выделена оранжевой линией). Это соответствует двум циклам, из которых собственно состоит исходный код обратного хода метода Гаусса. Из самого рисунка это заметить достаточно сложно, но анализ исходного кода показывает, что профиль образуется из обращений к 4 массивам. Три из них выделены на рис. 12.1 зеленым; к четвертому относятся все остальные обращения. Сразу отметим, что в данном профиле общее число обращений всего в несколько раз больше числа задействованных элементов массивов (4500 против 1000), а в таких условиях сложно добиться высокой локальности.<br />
<br />
Чтобы выяснить это, перейдем к подробному рассмотрению каждого из массивов в отдельности. <br />
<br />
На рис. 12.2 представлен фрагмент 1. Оранжевой линей также проведено разделение 2 этапов. Видно, что второй этап устроен очень просто и представляет собой последовательный перебор в обратном порядке. Первый этап имеет итерационный характер, причем на каждой итерации отбрасывается элемент массива с наименьшим номером. Подобный профиль характеризуется высокой пространственной локальностью, поскольку элементы перебираются подряд; достаточно высокой временной локальностью, поскольку на каждой итерации происходит повторное обращение к тем же элементам. <br />
<br />
[[file:gauss back 2.PNG|thumb|center|700px|Рисунок 12.2. Фрагмент 1 (профиль обращений к первому массиву)]]<br />
<br />
На рис. 12.3 представлен фрагмент 2, показывающий обращения ко второму массиву. Сразу заметим, что данный фрагмент еще меньше предыдущего – здесь задействовано всего 600 обращений в память. По сути, данный фрагмент является аналогичным этапу 1 фрагмент 1, с единственной разницей – здесь итерации расположены в обратном порядке. Однако это изменение не оказывает особого влияния ни на пространственную, ни не временную локальность, поэтому данный фрагмент обладает теми же свойствами.<br />
<br />
[[file:gauss back 3.PNG|thumb|center|700px|Рисунок 12.3. Фрагмент 2 (профиль обращений ко второму массиву)]]<br />
<br />
Далее рассмотрим фрагмент 3 (рис. 12.4). Здесь также оранжевой линией проведено разделение между двумя этапами, и также на второй этап приходится совсем немного обращений. Обращения на первом этапе образуют аналог случайного доступа, достаточно часто встречающийся, например, в случае косвенной адресации. При этом в некоторый момент к определенному элементу происходит достаточно много обращений подряд, после чего этот элемент более не используется. Такой профиль обычно характеризуется низкой пространственной и временной локальностью, что, однако, в данном случае в достаточной степени нивелируется малым числом задействованных элементов.<br />
<br />
[[file:gauss back 4.PNG|thumb|center|700px|Рисунок 12.4. Фрагмент 3 (профиль обращений к третьему массиву)]]<br />
<br />
Далее перейдем к фрагменту, занимающему основную часть рис. 12.1. Данный профиль отображен на рис. 12.5. В первую очередь стоит отметить особенность данного профиля – здесь задействовано более 1000 элементов, в то время как в остальных профилям – порядка 30. При этом число обращений даже меньше, чем обращений к первому или третьему массиву.<br />
<br />
Отметим также, что здесь большая часть обращений сосредоточена на втором этапе. Первый же этап является подобием первого этапа предыдущего массива (рис. 12.4) , за тем лишь исключением, что в данном профиле отсутствуют множественные обращения подряд к одному элементу перед тем, как элемент перестанет использоваться. С учетом того, что первый этап состоит всего из порядка 500 обращений, достаточно равномерно распределенных на отрезке в 1000 элементов, это говорит об очень низкой как пространственной, так и временной локальности.<br />
<br />
Другой характер обращений можно наблюдать на втором этапе. Видно, что он обладает более высокой пространственной локальностью, так как обращения явно сгруппированы в кластеры, при этом кластеры обладают подобным строением. Однако структура самого кластера видно плохо, поэтому требуется дальнейшее приближение.<br />
<br />
[[file:gauss back 5.PNG|thumb|center|700px|Рисунок 12.5. Фрагмент 4 (профиль обращений к четвертому массиву)]]<br />
<br />
На рис. 12.6 показаны два кластера, выделенные на рис. 12.5 зеленым цветом. Такой масштаб позволяет сразу увидеть, что каждый кластер представляет собой последовательный перебор с небольшим шагом определенного набора элементов массива.<br />
<br />
Следовательно, можно говорить о том, что второй этап фрагмента 4 обладает достаточно высокой пространственной локальностью (поскольку каждый кластер содержит последовательный перебор), но низкой временной локальностью (повторные обращения практически отсутствуют).<br />
<br />
[[file:gauss back 6.PNG|thumb|center|500px|Рисунок 12.6. Два кластера из фрагмента 4]]<br />
<br />
В целом по всему профилю можно сделать следующий вывод: первый три рассмотренных массива (особенно №1 и №2) обладают достаточно высокой локальностью, однако достаточно низкая пространственная и очень низкая временная локальность последнего массива в значительной степени снижают общую локальность всей программы.<br />
<br />
===== Количественная оценка локальности =====<br />
<br />
Основной фрагмент реализации, на основе которого были получены количественные оценки, приведен [http://git.algowiki-project.org/Voevodin/locality/blob/master/benchmarks/gauss_back/gauss_back.h здесь] (функция Kernel2). Условия запуска описаны [http://git.algowiki-project.org/Voevodin/locality/blob/master/README.md здесь].<br />
<br />
Первая оценка выполняется на основе характеристики daps, которая оценивает число выполненных обращений (чтений и записей) в память в секунду. Данная характеристика является аналогом оценки flops применительно к работе с памятью и является в большей степени оценкой производительности взаимодействия с памятью, чем оценкой локальности. Однако она служит хорошим источником информации, в том числе для сравнения с результатами по следующей характеристике cvg.<br />
<br />
На рисунке 12.7 приведены значения daps для реализаций распространенных алгоритмов, отсортированные по возрастанию (чем больше daps, тем в общем случае выше производительность). Можно увидеть, что, в соответствии с высказанным в предыдущем разделе предположением о негативном влиянии одного из массивов, данная программа показывает достаточно низкую производительность, заметно ниже, чем у прямого хода метода Гаусса.<br />
<br />
[[file:gauss back daps ru.PNG|thumb|center|700px|Рисунок 12.7. Сравнение значений оценки daps]]<br />
<br />
Вторая характеристика – cvg – предназначена для получения более машинно-независимой оценки локальности. Она определяет, насколько часто в программе необходимо подтягивать данные в кэш-память. Соответственно, чем меньше значение cvg, тем реже это нужно делать, тем лучше локальность. <br />
<br />
На рисунке 12.8 приведены значения cvg для того же набора реализаций, отсортированные по убыванию (чем меньше cvg, тем в общем случае выше локальность). Можно увидеть, что, согласно данной оценке, профиль обращений обладает низкой локальностью, лишь немногим лучше профиля программы со случайным доступом в память. Это повторяет выводы, сделанные на основе оценки daps.<br />
<br />
[[file:gauss back cvg.PNG|thumb|center|700px|Рисунок 12.8. Сравнение значений оценки cvg]]<br />
<br />
=== Возможные способы и особенности реализации параллельного алгоритма ===<br />
<br />
Вариантов параллельной реализации алгоритма не так уж и много, если не использовать то, что оба главных цикла можно развернуть, перейдя, таким образом, к блочной версии. Версии без развёртывания циклов возможны как с полностью параллельными циклами по I:<br />
<br />
<source lang="fortran"><br />
DO PARALLEL I = 1, N<br />
X(I) = Y(I)<br />
END DO<br />
DO J = N, 1, -1<br />
X(J) = X(J) / U(J,J)<br />
DO PARALLEL I = 1, J-1<br />
X(I) = X(I) - U(I,J)*X(J)<br />
END DO<br />
END DO<br />
</source><br />
<br />
так и с использованием "скошенного параллелизма" в главном гнезде циклов.<br />
<br />
=== Масштабируемость алгоритма и его реализации ===<br />
<br />
==== Описание масштабируемости алгоритма ====<br />
<br />
==== Описание масштабируемости реализации алгоритма ====<br />
<br />
=== Динамические характеристики и эффективность реализации алгоритма ===<br />
<br />
=== Выводы для классов архитектур ===<br />
<br />
Если исходить из структуры алгоритма, то при реализации на суперкомпьютерах следует выполнить две вещи. Во-первых, для минимизации обменов между узлами следует избрать блочный вариант, в котором или все элементы матрицы доступны на всех узлах, или заранее распределены по узлам. В таком случае количество передаваемых между узлами данных будет невелико по сравнению с количеством арифметических операций. Но при такой организации работы получится, что наибольшие временные затраты будут связаны с неоптимальностью обработки отдельных блоков. Поэтому, видимо, следует сначала оптимизировать не блочный алгоритм в целом, а подпрограммы, используемые на отдельных процессорах: точечный метод обратной подстановки, перемножения матриц и др. подпрограммы. Ниже содержится информация о возможном направлении такой оптимизации.<br />
<br />
=== Существующие реализации алгоритма ===<br />
<br />
Вещественный вариант обратной подстановки реализован как в основных библиотеках программ отечественных организаций, так и в западных пакетах LINPACK, LAPACK, SCALAPACK и др. <br />
При этом в отечественных реализациях, как правило, выполнены стандартные требования к методу с точки зрения ошибок округления, то есть, реализован режим накопления, и обычно нет лишних операций. <br />
Реализация точечного метода Холецкого в современных западных пакетах обычно происходит из одной и той же реализации метода в LINPACK, а та использует пакет BLAS. <br />
<br />
Для большинства пакетов существует блочный вариант обратной подстановки, в том числе и тот, граф которого топологически тождествен графу точечного варианта. Из-за того, что количество читаемых данных примерно равно количеству операций, блочность может дать некоторое ускорение работы благодаря лучшему использованию кэшей процессоров. Именно в направлении оптимизации кэширования и следует сосредоточить основные усилия при оптимизации работы программы.<br />
<br />
== Литература ==<br />
<br />
# В.В.Воеводин, Ю.А.Кузнецов. Матрицы и вычисления. - М.: Наука, 1984.<br />
# Дж.Форсайт, К.Моллер. Численное решение систем линейных алгебраических уравнений. - М.:Мир, 1969.<br />
<br />
<br />
[[En:Back substitution]]<br />
<br />
[[Категория:Статьи в работе]]</div>Oleg Arushanyanhttps://algowiki-project.org/w/ru/index.php?title=%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%A5%D0%BE%D0%BB%D0%B5%D1%86%D0%BA%D0%BE%D0%B3%D0%BE_(%D0%BD%D0%B0%D1%85%D0%BE%D0%B6%D0%B4%D0%B5%D0%BD%D0%B8%D0%B5_%D1%81%D0%B8%D0%BC%D0%BC%D0%B5%D1%82%D1%80%D0%B8%D1%87%D0%BD%D0%BE%D0%B3%D0%BE_%D1%82%D1%80%D0%B5%D1%83%D0%B3%D0%BE%D0%BB%D1%8C%D0%BD%D0%BE%D0%B3%D0%BE_%D1%80%D0%B0%D0%B7%D0%BB%D0%BE%D0%B6%D0%B5%D0%BD%D0%B8%D1%8F)&diff=1606Метод Холецкого (нахождение симметричного треугольного разложения)2015-04-16T13:32:17Z<p>Oleg Arushanyan: /* Точечный вариант */</p>
<hr />
<div>Основные авторы описания: [[Участник:Konshin|И.Н.Коньшин]]<br />
<br />
== Разложение Холецкого (метод квадратного корня), базовый точечный вещественный вариант для плотной симметричной положительно определённой матрицы ==<br />
<br />
=== <math>LL^T</math>-разложение ===<br />
<br />
Разложение Холецкого — представление симметричной положительно определённой матрицы <math>A=A^T>0</math> в виде произведения <math>A = LL^T</math>, где <math>L</math> — нижняя (''Lower'') треугольная матрица со строго положительными элементами на диагонали. Иногда разложение удобно записать в эквивалентной форме <math>A = U^TU</math>, где <math>U = L^T</math> — верхняя (''Upper'') треугольная матрица. Для любой симметричной положительно определённой матрицы разложение Холецкого существует и единственно.<br />
<br />
Элементы матрицы <math>L</math> можно вычислить, начиная с верхнего левого угла матрицы <math>A</math>, по формулам:<br />
<br />
:<math><br />
\begin{align}<br />
\ell_{ii} & = \sqrt{a_{ii} - \sum_{k=1}^{i-1} \ell_{ik}^2}, \\<br />
\ell_{ij} & = \frac{1}{\ell_{jj}} \left(a_{ij} - \sum_{k=1}^{j-1} \ell_{ik} \ell_{jk} \right), \quad j < i.<br />
\end{align}<br />
</math><br />
<br />
Выражение под квадратным корнем всегда положительно, если <math>A</math> — вещественная симметричная положительно определённая матрица.<br />
<br />
Вычисление происходит сверху вниз, слева направо, т.е. сначала вычисляется <math>L_{ij}</math> (<math>j < i</math>), а уже затем <math>L_{ii}</math>. Вычисления обычно проводятся в одной из следующих последовательностей.<br />
<br />
* Алгоритм Холецкого-Банашевича (''Cholesky–Banachiewicz algorithm'') или просто алгоритм Холецкого, когда вычисления начинаются с верхнего левого угла матрицы <math>L</math> и проводятся по строкам. Этот вариант разложения используется наиболее часто, особенно при использовании построчного формата хранения элементов матрицы <math>L</math>.<br />
<br />
* Краут-вариант алгоритма Холецкого (''Cholesky–Crout algorithm''), когда вычисления также начинаются с верхнего левого угла матрицы <math>L</math>, но проводятся по столбцам. Этот вариант разложения используется несколько реже, применяется он при использовании столбцевого формата хранения элементов матрицы <math>L</math>, а также когда необходимо проводить коррекцию ведущих элементов при выполнении приближенного разложения.<br />
<br />
Оба варианта разложения могут быть применены если требуется построить нижнетреугольный сомножитель <math>L</math> прямо поверх исходной матрицы <math>A</math>.<br />
<br />
В разделе [[Разложение Холецкого (метод квадратного корня)]] подробно рассмотрен базовый точечный вещественный вариант для плотной симметричной положительно определённой матрицы.<br />
<br />
=== <math>LDL^T</math>-разложение ===<br />
<br />
Иногда удобнее бывает рассматривать <math>LDL^T</math> вариант симметричного треугольного разложения, в котором матрица <math>L</math> является нижней унитреугольной (т.е. имеет единицы на главной диагонали), а <math>D</math> - диагональная матрица с положительными элементами. В этом варианте разложения легко проследить связь как с ранее рассмотренным <math>LL^T</math> вариантом:<br />
<br />
:<math>A = LDL^T = LD^{1/2}D^{1/2}L^T = (LD^{1/2})\,(LD^{1/2})^T = \tilde L \tilde L^T,</math><br />
<br />
так и с несимметричным <math>LU</math>-разложением:<br />
<br />
:<math>A = LDL^T = L(DL^T) = LU.</math><br />
<br />
== Разложение Холецкого, блочный вещественный вариант для плотной симметричной положительно определённой матрицы ==<br />
<br />
Можно также рассмотреть блочный вариант разложения Холецкого. Предположим, что <math>n=MN</math>, тогда исходную матрицу <math>A</math> размера <math>n\times n</math> можно представить как блочную матрицу размера <math>N\times N</math> с блоками размера <math>M\times M</math>.<br />
Все формулы, используемые для получения точечного разложения Холецкого, для блочной матрицы <math>A</math> останутся практически без изменений.<br />
Вместо явного обращения диагональных блоков, эффективнее хранить их в факторизованном виде <math>D_{ii}=L_{ii}L^T_{ii}</math>,<br />
а вместо операции деления использовать соответствующие операции решения для треугольных систем.<br />
Общее количество арифметических операций при этом останется практически неизменным, но зато существенно возрастет локальность вычислений.<br />
Размер блока <math>M</math> выбирают таким образом, чтобы все блоки, участвующие в операции исключения, помещались в кэш первого или второго уровня. В этом случае подкачки данных в память будут минимальными.<br />
<br />
Аналогичный прием понадобится также и для эффективной реализации параллельной версии разложения Холецкого, что позволит минимизировать как общее количество межпроцессорных обменов, так и количество пересылаемой между процессорами информации.<br />
Полезным побочным эффектом применения блочной версии разложения Холецкого может стать повышение скалярной эффективности алгоритма за счет явного использования размера блока <math>M</math> во внутренних циклах (прием "разворачивание цикла" или "loop unrolling").<br />
<br />
== Разложение Холецкого, точечный вещественный вариант для разреженной симметричной положительно определённой матрицы ==<br />
<br />
Если исходная матрица <math>A</math> представлена в разреженном виде, то для экономии памяти, а также арифметических операций, необходимо учитывать ее разреженность.<br />
<br />
=== Основные отличия от случая плотной матрицы ===<br />
<br />
В этом разделе необходимо рассмотреть матрицы, характеризующиеся способом хранения ненулевых элементов, и имеющие следующие виды разреженности.<br />
<br />
* ''Лента'' - матрица, ненулевые элементы которой сосредоточены внутри ленты шириной <math>2d+1</math>, т.е. когда <math>a_{ij}=0</math> при <math>|i-j| > d</math>. В этом случае, при проведении разложения Холецкого новые ненулевые элементы могут образовываться только внутри этой же ленты. Количество ненулевых элементов в исходной матрице <math>A</math>, а также в нижнетреугольном множителе <math>L</math> будет около <math>(d+1)n</math>, а арифметические затраты составят приблизительно <math>d^2n</math>.<br />
<br />
* ''Профиль'' - в более общем случае, заполнение в каждой строке треугольного множителе <math>L</math> будет определяться позицией первого ненулевого элемента. Сумма по всем строкам расстояний от первого ненулевого элемента строки до главной диагонали и составляет "профиль" матрицы и определяет верхнюю границу количества ненулевых элементов в нижнетреугольном множителе <math>L</math>.<br />
<br />
* ''Общая структура разреженности''. Верхней границей заполнения треугольного множителя <math>L</math>, конечно же, будет значение "профиля" матрицы, но учет особенностей структуры ненулевых элементов внутри профиля иногда может дать дополнительный эффект в повышении эффективности вычислений.<br />
<br />
При рассмотрении общего случая разреженности необходимо выбрать формат хранения разреженных данных. Таковым может быть, например, формат построчного сжатия данных ("compressed sparse row" или CSR формат). В первом вещественном массиве, подряд (обычно в порядке возрастания номеров столбцов) хранятся ненулевые элементы матрицы, во втором, в том же порядке хранятся номера столбцов, в третьем, отдельно сохраняется начало каждой строки. Если общее количество ненулевых элементов в матрице равно nnz ("number of nonzeros"), то память для хранения разреженных данных такой матрицы в формате CSR при использовании двойной точности составит <math>3\,{\rm nnz}+n+1</math>. Оценку количества арифметических операций в общем случае невозможно, т.к. помимо количества ненулевых элементов в исходной матрице оно существенно зависит от структуры ее разреженности.<br />
<br />
Для реализации разложения Холецкого в этом случае понадобится несколько операций с разреженными строками:<br />
<br />
* ''копирование'' из одной разреженной строки в другую (или во временный "плотный" вектор, операция ''распаковки'' данных);<br />
<br />
* выполнение ''операции исключения'' для одного из элементов строки;<br />
<br />
* ''вставка'' в строку нового ненулевого элемента ("fill-in");<br />
<br />
* ''сжатие'' данных с копированием из временного плотного вектора в сжатый разреженный (операция ''упаковки'' данных).<br />
<br />
=== Переупорядочивания для уменьшения количества новых ненулевых элементов ===<br />
<br />
Структура треугольного множителя <math>L</math>, а также объем памяти им занимаемый зависят от упорядочивания строк и столбцов исходной матрицы <math>A</math>, в котором проводилось разложение. Существуют алгоритмы, минимизирующие заполнение матрицы <math>L</math>.<br />
<br />
* В первую очередь это алгоритм RCM (Reverse Cuthill–McKee), который предназначен для уменьшения профиля матрицы. Одновременно с уменьшением профиля происходит и уменьшение заполнения треугольного множителя <math>L</math>. Это очень широко применяемый, быстрый, но не самый эффективный алгоритм.<br />
<br />
* Алгоритм вложенных сечений (Nested Dissection, ND) - служит именно для минимизации заполнения множителя <math>L</math>. В некоторых частных случаях доказана его асимптотическая оптимальность.<br />
<br />
В общем случае, проблема поиска перестановки, минимизирующей заполнение множителя <math>L</math>, является NP-полной задачей.<br />
<br />
== Разложение Холецкого, блочный вещественный вариант для разреженной симметричной положительно определённой матрицы ==<br />
<br />
Иногда разреженную симметричную матрицу бывает удобно представить в блочном виде с блоками небольшого размера <math>M</math>, равного, например, количеству неизвестных функций на узел при конечно-элементной или конечно-разностной аппроксимации уравнений в частных производных.<br />
В этом случае структура разреженности хранится для всей блочной структуры разреженности (что позволяет экономить память на хранении целочисленных массивов).<br />
Если общее количество ненулевых блоков размера <math>M\times M</math> в матрице равно nnz ("number of nonzeros"), то память для хранения разреженных данных такой мелкоблочной матрицы в формате CSR при использовании двойной точности составит <math>(2M^2+1)\,{\rm nnz}+n/M+1</math>.<br />
<br />
В некоторых случаях, размер блока <math>M</math> может выбираться из других соображений, например, для повышения эффективности работы процедур нижнего уровня за счет приема разворачивания циклов (''loop unrolling'').<br />
<br />
Алгоритмы, необходимые при выполнении разложения Холецкого для матриц, рассмотренных в этом разделе, могут быть получены комбинацией уже рассмотренных идей <br />
[[Метод_Холецкого_(нахождение_симметричного_треугольного_разложения)#Разложение Холецкого, блочный вещественный вариант для плотной симметричной положительно определённой матрицы|блочности]] и [[Метод_Холецкого_(нахождение_симметричного_треугольного_разложения)#Разложение Холецкого, точечный вещественный вариант для разреженной симметричной положительно определённой матрицы|разреженности]].<br />
<br />
== Разложение Холецкого для симметричной незнакоопределенной (седловой) матрицы ==<br />
<br />
Если симметричная матрица <math>{\mathcal A}</math> представима в виде<br />
<br />
:<math><br />
{\mathcal A} =<br />
\begin{bmatrix}<br />
A & B^T \\ <br />
B & -C <br />
\end{bmatrix} ,<br />
</math><br />
<br />
где <math>A</math> - симметричная положительно определенная (<math>A=A^T>0</math>) и <math>C</math> - симметричная неотрицательно определенная (<math>C=C^T\ge0</math>) матрицы, то, выполнив один шаг блочного исключения, ее можно преобразовать к виду<br />
<br />
:<math><br />
\begin{bmatrix}<br />
A & 0 \\ <br />
0 & S <br />
\end{bmatrix} ,<br />
</math><br />
<br />
где матрица дополнения по Шуру <math>S=-(C+B^TA^{-1}B)</math> является строго отрицательно определенной (<math>S=S^T<0</math>).<br />
Это означает, что матрица <math>{\mathcal A}</math> имеет <math>n_A</math> положительных и <math>n_C</math> отрицательных собственных значений, где через <math>n_A</math> и <math>n_C</math> обозначены размерности матриц <math>A</math> и <math>C</math>, соответственно.<br />
<br />
В этом случае существует симметричное треугольное разложение вида <math>{\mathcal A}={\mathcal L}D{\mathcal L}^T</math>, где <math>{\mathcal L}</math> является нижней унитреугольной, а диагональная матрица <math>D</math> содержит <math>n_A</math> положительных и <math>n_C</math> отрицательных элементов на главной диагонали, причем такое разложение может быть получено напрямую без выбора ведущего элемента даже если <math>C</math> - нулевая матрица.<br />
<br />
В общем случае разложения невырожденной незнакоопределенной системы необходимо применять выбор ведущего элемента с главной диагонали матрицы, что соответствует некоторой симметричной перестановке строк и столбцов исходной матрицы <math>{\mathcal A}</math>.<br />
<br />
== Разложение Холецкого для эрмитовой матрицы ==<br />
<br />
Эрмитовой (или комплексно-самосопряженной) матрицей называют такую квадратную комплексную матрицу <math>A</math>, для элементов которой выполняется соотношение<br />
<math>a_{ij}=\overline{a_{ji}}</math> (здесь, если <math>z=a+{\rm i\,}b\,</math> и <math>{\rm i}^2=-1</math>, то <math>\overline z=a-{\rm i\,}b\,</math>). В матричном виде это можно записать как <math>A=\overline{A^T}</math> или <math>A=A^*=A^Н</math>.<br />
<br />
=== Точечный вариант ===<br />
<br />
Как естественное обобщение точечного разложения Холецкого для симметричной положительно определеной матрицы может быть рассмотрено разложение Холецкого для эрмитовой положительно определеной матрицы. Все формулы для вычисления разложения остаются прежними, только теперь вместо операций над вещественными числами выполняются аналогичные комплексные операции:<br />
<br />
:<math><br />
\begin{align}<br />
L_{ii} & = \sqrt{ A_{ii} - \sum_{k=1}^{i-1} L_{ik}L_{ik}^* }, \\<br />
L_{ij} & = \frac{1}{L_{jj}} \left( A_{ij} - \sum_{k=1}^{j-1} L_{ik} L_{jk}^* \right), \quad j < i.<br />
\end{align}<br />
</math><br />
<br />
В отличие от<br />
[[Метод_Холецкого_(нахождение_симметричного_треугольного_разложения)#Разложение Холецкого (метод квадратного корня), базовый точечный вещественный вариант для плотной симметричной положительно определённой матрицы|вещественного варианта]],<br />
при выполнении аналогичных комплексных операций потребуется считывать из памяти вдвое больше данных и производить над ними примерно вчетверо больше арифметических операций, что должно не только несколько улучшить локальность вычислений, но и повысить общую эффективность вычислений.<br />
<br />
=== Блочный вариант ===<br />
<br />
Реализация блочного варианта разложения Холецкого для эрмитовых матриц аналогична рассмотренному выше<br />
[[Метод_Холецкого_(нахождение_симметричного_треугольного_разложения)#Разложение Холецкого, блочный вещественный вариант для плотной симметричной положительно-определённой матрицы|блочному варианту]]<br />
для вещественных матриц.<br />
<br />
== Использование разложения Холецкого в итерационных методах ==<br />
<br />
При выполнении разложения Холецкого в арифметике с фиксированной машинной точностью полученные треугольный фактор <math>L</math> и само решение может оказаться недостаточно точным.<br />
Для получения более точного решения может применяться некоторый итерационный метод (например, метод сопряженных градиентов), с использованием полученного разложения <math>LL^T</math> в качестве предобусловливателя.<br />
<br />
Основной причиной формирование неполного или неточного разложения в качестве предобусловливателя чаще всего бывает требование экономии памяти.<br />
<br />
=== Ограничивание заполнения в разложении Холецкого ===<br />
<br />
При выполнении разложения Холецкого для разреженной матрицы, может образовываться такое большое количество новых ненулевых элементов, что оперативной памяти на хранение полного разложения окажется недостаточно.<br />
В этом случае можно построить неполное или приближенное разложение для применения его в дальнейшем в качестве предобусловливателя.<br />
В англоязычной литературе для обозначения таких разложения применяют единый термин Incomplete Cholesky factorization, или сокращенно IC-разложение.<br />
<br />
=== Неполное разложение Холецкого по позициям IC(<math>k</math>) ===<br />
<br />
Неполное разложение Холецкого можно получить используя заранее выбранные ограничения по структуре заполнения.<br />
Чаще всего получают разложение Холецкого на тех же позициях, в которых находятся ненулевые элементы исходной матрицы <math>A</math>. Такое разложение обозначают IC(0) или просто IC0.<br />
<br />
Если качества разложения IC0 оказывается недостаточно, то можно выбрать более широкую структуру треугольного множителя <math>L</math>, например, разрешить образование одного уровня новых ненулевых элементов от исходной структуры матрицы <math>A</math>. Формально, это означает заполнение внутри структуры матрицы <math>A^2</math>, а такое разложение обозначают IC(1).<br />
<br />
Можно рассмотреть и более общий случай, с заполнением внутри структуры матрицы <math>A^{k+1}</math>, где <math>k \geq 0</math>. Такое разложение обозначают IC(<math>k</math>).<br />
<br />
Обычно с ростом значения <math>k</math> точность неполного разложения IC(<math>k</math>) возрастает, хотя это совсем не является обязательным даже для симметричных положительно определенных матриц, полное разложение для которых существует и находится однозначно.<br />
Из-за неполноты разложения на главной диагонали могут оказаться нулевые или даже отрицательные элементы.<br />
Чтобы завершить треугольное разложение в таких случаях применяют предварительный диагональный сдвиг исходной матрицы <math>A+\varepsilon I</math> перед ее разложением. Здесь <math>\varepsilon>0</math> - малый параметр, а <math>I</math> - диагональная матрица.<br />
Если слишком малый или не положительный диагональный элемент образуется в процессе разложения, то применяют его замену на некоторое заранее выбранное значение. Такую операцию называют диагональной коррекцией разложения.<br />
<br />
Неполное разложение IC(<math>k</math>) иногда называют также "разложение по позициям".<br />
<br />
=== Приближенное разложение Холецкого по значениям IC(<math>\tau</math>) ===<br />
<br />
Для контроля заполнения в треугольном множителе <math>L</math> разложения Холецкого, кроме структурных ограничений, можно также применить ограничение разложения в зависимости от значения самих элементов разложения.<br />
Например, можно сохранять только элементы, большие по модулю чем некоторый малый параметр <math>\tau>0</math>.<br />
В этом случае разложение называют ''приближенным'' разложением Холецкого или разложением "по значению" и обозначают IC(<math>\tau</math>).<br />
Величину <math>\tau</math> называют "порогом" разложения или "порогом" фильтрации.<br />
<br />
Вполне правомерным является ожидание того, что в уменьшением <math>\tau</math> точность полученного разложения будет возрастать, правда за счет роста количества ненулевых элементов в треугольном множителе <math>L</math>.<br />
Недостатком же такого разложения является то, что, в общем случае, предсказать заполнение <math>L</math> не возможно.<br />
<br />
С точки зрения устойчивости разложения вариант приближенного разложения Холецкого является более предпочтительным, хотя применение предварительного диагонального сдвига, а также диагональной коррекции также допускается.<br />
Если же описанные приемы не помогаю получить разложения достаточной точности, то можно применить прием модификации диагонали Азиза-Дженингса, который при отбрасывании малого элемента разложения <math>\ell_{ij}</math> состоит в добавлении модуля этого элемента к диагональным элементам разложения <math>\ell_{ii}</math> и <math>\ell_{jj}</math>. Это прием гарантирует существование приближенного разложения для любой симметричной положительно определенной матрицы <math>A</math>. Наиболее эффективно этот прием модификации главной диагонали можно организовать при использовании Ктаут-версии разложения Холецкого.<br />
<br />
=== Приближенное разложение Холецкого второго порядка IC(<math>\tau_1,\tau_2</math>) ===<br />
<br />
Для повышения точности приближенного разложения можно применить "двухпороговую" версию приближенного разложения Холецкого. Основная идея такого разложения, называемого разложением Тисменецкого-Капорина, состоит в том чтобы вычисление разложения проводить в более высокой точности <math>\tau_2</math>, а сохранять в треугольном множителе только значения, которые по модулю не меньше <math>\tau_1</math>. Обычно полагают <math>\tau_1=\tau</math> и <math>\tau_2=\tau^2</math>, в этом случае разложение называют разложением "второго порядка", т.к. элементы матрицы ошибок оказываются по модулю меньше чем <math>\tau^2</math>.<br />
Для обозначения симметричного порогового разложения второго порядка используют обозначение IC2, которое не следует путать со структурным разложением IC(2) (т.е. разложением IC(<math>k</math>), где <math>k=2</math>).<br />
<br />
Такое разложение обычно используется вместе с приемом Азиза-Дженингса для модификации диагональных элементов, получая вариант "безотказного" разложения для любой симметричной положительно определенной матрицы <math>A</math>.<br />
Этот вариант разложения позволяет получать наиболее точные разложения (при одинаковом заполнении множителя <math>L</math>), хотя для их вычисления приходится тратить больше времени на вычисление самого разложения.<br />
<br />
=== Комбинация разложений Холецкого IC(<math>k,\tau</math>) и IC(<math>\tau,m</math>) ===<br />
<br />
Для экономии памяти при вычислении неполного или приближенного разложения Холецкого можно использовать следующие два варианта симметричных треугольных разложений.<br />
<br />
Для контроля верхней границы заполнения треугольного множителя <math>L</math> можно предложить использовать заполнение как и для разложения IC(<math>k</math>), при некотором выбранном значении <math>k</math>. Для дальнейшей экономии памяти разложение в заданной структуре разреженности можно вести с использованием порога разложения <math>\tau</math>, как и при проведении разложения IC(<math>\tau</math>).<br />
Такую комбинацию можно назвать IC(<math>k,\tau</math>)-разложением.<br />
Применяться она может, например, при необходимости априорных структурных ограничений для минимизации обменов при использовании параллельной версии разложения для распределенной памяти.<br />
<br />
Второй вариант структурно-порогового разложения можно описать следующим образом.<br />
При проведении обычного порогового IC(<math>\tau</math>) разложения, наложим дополнительное ограничение на элементы строк матрицы <math>L</math>: разрешим сохранение только не более чем <math>m</math> наибольших по модулю элементов рассматриваемой строки множителя <math>L</math>.<br />
При общей размерности задачи <math>n</math> в матрице <math>L</math> будет не более чем <math>nm</math> элементов.<br />
Такой подход представляется разумным, например, для матриц полученных в результате дискретизации с достаточно регулярным шаблоном.<br />
Наиболее известен несимметричный вариант такого разложения, предложенного Саадом и называемого ILUT-разложением.<br />
<br />
== Использование разложения Холецкого в параллельных итерационных алгоритмах ==<br />
<br />
Формулы разложения Холецкого по большей части имеют рекурсивный характер и выделение параллельных и независимых этапов вычислений является не очевидной и непростой задачей.<br />
Слишком прямолинейное ее решение может привести к значительному объему пересылаемых данных, что значительно снизит результат распараллеливания.<br />
Наибольший эффект может дать подход, основанный на предварительном переупорядочивании исходной матрицы.<br />
<br />
=== Переупорядочивания для выделения блочности ===<br />
<br />
Для того чтобы выделить независимые блоки вычислений, можно использовать симметричные перестановки строк и столбцов исходной матрицы, приводящие ее к блочно окаймленному виду.<br />
В этом случае основная часть работы будет сосредоточена в независимых блоках, которые могут обрабатываться параллельно и без обменов данными.<br />
<br />
Наиболее простым, но не слишком эффективных способом упорядочивания является предварительное упорядочивание матрицы с помощью обратного алгоритма Катхилла-Макки (RCM, reverse Cuthill—McKee) для минимизации ширины профиля, а затем равномерное разбиение по блокам (процессорам) в новом упорядочивании. После присваивания номера процессора каждой вершине графа матрицы, в независимые блоки можно выделить те вершины графа, которые связаны только с вершинами, имеющими тот же номер процессора (т.е. являющимися внутренними вершинами). Остальные вершины можно объединить в последний блок окаймления, который будет обрабатываться отдельно. Все вычисления внутри блоков будут независимы и могут выполняться параллельно.<br />
Для повышения эффективности треугольной факторизации внутренние вершины каждого из блоком можно также упорядочить с помощью метода RCM.<br />
<br />
Более эффективными с точки зрения минимизации ширины окаймления будут следущие методы:<br />
<br />
* Метод минимальных сепараторов, который заключается в последовательном нахождении имеющих минимальный размер разделителей (сепараторов), обеспечивающих расщепление оставшихся вершин на два независимых блока.<br />
<br />
* Метод минимальной степени (Minimum Degree, MD). Прямое применение этого метода к матрицах большого размера затруднительно, поэтому используется ''приближенный метод минимальной степени'' (Approximate Minimum Degree, AMD).<br />
<br />
* Метод вложенных сечений (Nested Dissection, ND). Это рекурсивный алгоритм, на каждом шаге разделяющий множество вершин на два независимых блока, представленые в <math>2\times2</math> блочно-окаймленном виде.<br />
В качестве побочного положительного эффекта от такого упорядочивания, для некоторого вида матриц доказано, что полное разложение будет иметь почти минимально возможное количество ненулевых элементов.<br />
Нахождение оптимальной перестановки в общем случае является NP-полной задачей.<br />
<br />
Существуют и другие алгоритмы упорядочивания матриц для наиболее оптимального их распределения по процессорам.<br />
Наиболее популярными являются последовательные пакеты METIS, JOSTLE, SCOTCH, CHACO, PARTY, а также параллельные коды PARMETIS, JOSTLE, PT-SCOTCH и ZOLTAN.<br />
Многие из них являются свободно распространяемыми.<br />
<br />
=== Разложение в независимых блоках ===<br />
<br />
Вычисления в независимых блоках полностью независимы и могут выполняться параллельно без обменов данными между процессорами.<br />
Единственным недостатком может быть лишь то, что для пороговых разложений количество арифметических операций на различных процессорах может различаться, что может привести к некоторой несбалансированности вычислений.<br />
<br />
=== Разложение в сепараторах ===<br />
<br />
Последний блок, в котором собраны сепараторы всех блоков при небольшом количестве используемых процессоров может обрабатываться, например, на одном головном процессоре.<br />
Если же процессоров достаточно много, но обработку сепараторов также необходимо производить совместно.<br />
<br />
=== Иерархические и вложенные алгоритмы ===<br />
<br />
Для обработки сепараторов могут быть применены те же алгоритмы упорядочивания и выделения независимых блоков что и для исходной задачи.<br />
Этот же прием может быть применен и на следующем шаге, с получением иерархического или вложенного алгоритма.<br />
<br />
В случае применения порогового разложения такие построения могут быть явно применены к построенному дополнению по Шуру.<br />
Для структурных факторизаций существуют приложения сразу строящие многоуровневые упорядочивания и обеспечивающие минимальность заполнения и обменов.<br />
<br />
=== Блочный метод Якоби ===<br />
<br />
Описанные выше подходы относились к типу "прямого" или "явного" разложения, когда при отбрасывании элементов разложения в расчет принимались только их абсолютные значения. Структурные свойства при этом являлись подчиненными и на отбрасывание элементов никак не влияли.<br />
<br />
Альтернативным является подход, когда из-за соображений увеличения ресурса параллелизма некоторые элементы отбрасываются исключительно из структурных соображений. Например, можно сначала распределить строки матрицы по процессорам, а затем перед проведением разложения просто отбросить все элементы связывающие один процессор с другим.<br />
В этом случае разложение будет проходить полностью независимо в каждом из блоков.<br />
Внутри каждого блока может проводиться любой из структурных или пороговых вариантов разложения Холецкого.<br />
Построение такого предобусловливателя называют блочным методом Якоби без перекрытия блоков (Block Jacobi, BJ).<br />
Такое предобусловливание является наиболее простым, применение его наиболее параллельным (полностью без обменов данными), правда сходимость может оставлять желать лучшего, почти независимо от качества разложения внутри каждого из блоков.<br />
<br />
=== Аддитивный метод Шварца ===<br />
<br />
Разложение гораздо более высокого качества по сравнению с методом Якоби можно получить применяя аддитивный метод Шварца (Additive Schwarz, AS). Иногда этот метод называют также методом Якоби с перекрытиями. Суть его заключается в расширении структуры каждого из блоков матрицы за счет добавления нескольких слоев близлежащих строк матрицы. Треугольное разложение строится для расширенной матрицы, таким образом на каждом из процессоров происходит решение расширенной подзадачи с привлечением данных от других процессоров.<br />
После нахождения решения подзадачи на каждом из процессоров обычно происходит отбрасывание не локальных компонент решения. Такой вариант метода называют аддитивный метод Шварца с ограничениями (Restricted Additive Schwarz, RAS).<br />
<br />
Сходимость аддитивного метода Шварца бывает гораздо выше чем сходимость метода Якоби, и обычно монотонно улучшается с ростом размера перекрытия. Несмотря на дополнительные вычисление и обмены, общее время решения на параллельном компьютере может быть существенно меньше.<br />
<br />
=== Неполное обратное треугольное разложения ===<br />
<br />
Существует и другой вариант аддитивного разложения, который кроме несколько более быстрой сходимости опирается на построение перекрытий блоков только в одну сторону ("назад", т.е. в сторону меньших номеров строк).<br />
Название этого метода блочное неполного обратного разложения Холецкого, имеющее только английскую аббревиатуру BIIC (Block Incomplete Inverse Cholesky).<br />
Позднее, вместе с рассмотрением несимметричного варианта разложения (BIILU), этот метод стал называться методом неполного обратного треугольного разложения, или НОТ-разложения.<br />
<br />
Комбинация этого метода с неполным симметричным треугольным разложением второго порядка IC2 внутри каждого из блоков имеет обозначение BIIC2.<br />
<br />
Идея этого метода впервые была предложена И.Е.Капориным в виде последовательного алгоритма.<br />
В литературе встречается также название этого метода как метод Капорина-Коньшина, по имени авторов впервые представивших его параллельную реализацию и проанализировавших ее свойства.<br />
<br />
== Решение линейных систем с треугольной матрицей ==<br />
<br />
Разложение Холецкого может применяться для решения системы линейных уравнений <math>Ax = b</math>, если матрица <math>A</math> симметрична и положительно определена. <br />
Выполнив разложение <math>A = LL^T</math>, решение <math>x</math> получается последовательным решением двух треугольных систем уравнений <math>Ly = b</math> и <math>L^T x = y</math>.<br />
<br />
=== Решение системы с плотной нижнетреугольной матрицей ===<br />
<br />
Решение линейной системы с плотной нижнетреугольной матрицей <math>L y = b</math> можно представить в виде "прямого" хода, т.е. цепочки вычислений, начиная с верхнего угла матрицы <math>L</math> по возрастанию номера строки <math>i</math>:<br />
<br />
:<math><br />
\begin{align}<br />
y_{1} & = b_{1}, \\<br />
y_{i} & = b_{i} - \sum_{j = 1}^{i-1} \ell_{ij} y_{j}, \quad i = 2,...,n.<br />
\end{align}<br />
</math><br />
<br />
В разделе<br />
[[Прямая_подстановка_(вещественный_вариант)]]<br />
содержится подробное описание алгоритма и его анализ.<br />
<br />
=== Решение системы с плотной верхнетреугольной матрицей ===<br />
<br />
Решение линейной системы с плотной верхнетреугольной матрицей <math>U x = y</math> (где, например, <math>U=L^T</math>)<br />
можно представить в виде "обратного" хода, т.е. цепочки вычислений, начиная с нижнего угла матрицы <math>U</math> при убываниии номера строки <math>i</math>:<br />
<br />
:<math><br />
\begin{align}<br />
x_{n} & = y_{n}/u_{nn}, \\<br />
x_{i} & = \left (y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} \right ) / u_{ii}, \quad i = n - 1,...,1.<br />
\end{align}<br />
</math><br />
<br />
В разделе<br />
[[Обратная_подстановка_(вещественный_вариант)]]<br />
содержится подробное описание алгоритма и его анализ.<br />
<br />
=== Решение системы с разреженной нижнетреугольной матрицей ===<br />
<br />
Решение линейных систем с разреженной нижне- или верхнетреугольной матрицей аналогично рассмотренным алгоритмам для плотных матриц, при этом подстановки ведутся исключительно для ненулевых элементов с учетом идеи работы с <br />
[[Метод_Холецкого_(нахождение_симметричного_треугольного_разложения)#Разложение Холецкого, точечный вещественный вариант для разреженной симметричной положительно определённой матрицы|разреженными]]<br />
матрицами.<br />
<br />
=== Решение системы с комплексной треугольной матрицей ===<br />
<br />
Решение линейных систем с комплексной нижне- или верхнетреугольной матрицей аналогично рассмотренным алгоритмам для вещественных матриц, при этом арифметические операции выполняются в комплексной арифметике, аналогично операциям раздела <br />
[[Метод_Холецкого_(нахождение_симметричного_треугольного_разложения)#Разложение Холецкого для эрмитовой матрицы|факторизации эрмитовых матриц]].<br />
<br />
=== Решение систем с блочноокаймленными треугольными матрицами ===<br />
<br />
Особенность решения линейных систем с блочноокаймленными треугольными матрицами в том что независимость вычислений в отдельных блоках дает возможность проведения параллельных вычислений.<br />
<br />
== Существующие реализации алгоритма ==<br />
<br />
* В [http://netlib.org/lapack/ LAPACK] используется функция [http://www.netlib.org/lapack/double/dpotrf.f DPBTRF] (последовательная реализация для двойной точности).<br />
<br />
* В [http://www.netlib.org/scalapack/ ScaLAPACK] используется функция [http://www.netlib.org/scalapack/double/pdpbtrf.f PDPBTRF] (параллельная реализация для двойной точности).<br />
<br />
* В [https://ru.wikipedia.org/wiki/SAS_Institute SAS] используется функция ROOT( matrix ), входящая в пакет SAS IML.<br />
<br />
* В системах [https://ru.wikipedia.org/wiki/MATLAB MATLAB], [https://ru.wikipedia.org/wiki/GNU_Octave Octave], [https://ru.wikipedia.org/wiki/R_(язык_программирования) R] разложение выполняется командой U = chol(A).<br />
<br />
* В [https://ru.wikipedia.org/wiki/Maple Maple] и [https://ru.wikipedia.org/wiki/NumPy NumPy] существует процедура cholesky в модуле linalg.<br />
<br />
* В [https://ru.wikipedia.org/wiki/Mathematica Mathematica] используется процедура CholeskyDecomposition[A].<br />
<br />
* В [https://ru.wikipedia.org/wiki/GSL GSL] используется функция gsl_linalg_cholesky_decomp.<br />
<br />
* В [http://www.alglib.net ALGLIB] имеются реализации как [http://www.alglib.net/matrixops/cholesky.php LLT] так и [http://www.alglib.net/matrixops/symmetric/ldlt.php LDLT] разложений для различных языков программирования: C#, C++, C++ (арифметика повышенной точности), FreePascal, Delphi, VB.NET, VBA, Python.<br />
<br />
* В [http://www.bluebit.gr/matrix-calculator/ Online Matrix Calculator] непосредственно в web-интерфейсе можно выполнить разложение Холецкого, выбрав раздел ''Cholesky Decomposition''.<br />
<br />
<br />
[[en:Cholesky method]]<br />
<br />
[[Категория:Статьи в работе]]</div>Oleg Arushanyan