Метод Холецкого (нахождение симметричного треугольного разложения)

Основные авторы описания: И.Н.Коньшин

1 Разложение Холецкого (метод квадратного корня), базовый точечный вещественный вариант для плотной симметричной положительно-определённой матрицы

1.1 [math]LL^T[/math] разложение

Разложение Холецкого — представление симметричной положительно-определённой матрицы [math]A[/math] в виде произведения [math]A = LL^T[/math], где [math]L[/math] — нижняя (Lower) треугольная матрица со строго положительными элементами на диагонали. Иногда разложение удобно записать в эквивалентной форме [math]A = U^TU[/math], где [math]U = L^T[/math] — верхняя (Upper) треугольная матрица. Разложение Холецкого всегда существует и единственно для любой симметричной положительно-определённой матрицы.

Элементы матрицы [math]L[/math] можно вычислить, начиная с верхнего левого угла матрицы, по формулам:

[math] \begin{align} \ell_{ii} & = \sqrt{a_{ii} - \sum_{k=1}^{i-1} \ell_{ik}^2}, \\ \ell_{ij} & = \frac{1}{\ell_{jj}} \left(a_{ij} - \sum_{k=1}^{j-1} \ell_{ik} \ell_{jk} \right), \quad j \lt i. \end{align} [/math]

Выражение под квадратным корнем всегда положительно, если [math]A[/math] — действительная положительно-определённая матрица.

Вычисление происходит сверху вниз, слева направо, т.е. сначала вычисляется [math]L_{ij}[/math] ([math]j \lt i[/math]), а уже затем [math]L_{ii}[/math]. Вычисления обычно проводятся в одной из следующих последовательностей:

• алгоритм Холецкого-Банашевича (Cholesky–Banachiewicz algorithm) или просто алгоритм Холецкого, когда вычисления начинаются с верхнего левого угла матрицы [math]L[/math] и проводятся по строкам. Этот вариант разложения используется наиболее часто, особенно при использовании построчного формата хранения элементов матрицы [math]L[/math].

• Краут-вариант алгоритма Холецкого (Cholesky–Crout algorithm), когда вычисления также начинаются с верхнего левого угла матрицы [math]L[/math], но проводятся по столбцам. Этот вариант разложения используется несколько реже, применяется он при использовании столбцевого формата хранения элементов матрицы [math]L[/math], а также когда необходимо проводить коррекцию ведущих элементов при выполнении приближенного разложения.

Оба варианта разложения могут быть применены если требуется построить нижнетреугольный сомножитель [math]L[/math] прямо поверх исходной матрицы [math]A[/math].

В разделе Разложение Холецкого (метод квадратного корня) подробно рассмотрен базовый точечный вещественный вариант для плотной симметричной положительно-определённой матрицы.

1.2 [math]LDL^T[/math] разложение

Иногда удобнее бывает рассматривать [math]LDL^T[/math] вариант симметричного треугольного разложения, в котором матрица [math]L[/math] является унитреугольной (т.е. имеет единицы на главной диагонали), а [math]D[/math] - диагональная матрица с положительными элементами. В этом варианте разложения легко проследить связь как с ранее рассмотренным [math]LL^T[/math] вариантом:

[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]

так и с несимметричным [math]LU[/math] разложением:

[math]A = LDL^T = L(DL^T) = LU.[/math]

2 Разложение Холецкого, блочный вещественный вариант для плотной симметричной положительно-определённой матрицы

Можно также рассмотреть блочный вариант разложения Холецкого. Предположим, что [math]n=MN[/math], тогда исходную матрицу [math]A[/math] размера [math]n\times n[/math] можно представить как блочную матрицу размера [math]N\times N[/math] с блоками размера [math]M\times M[/math]. Все формулы, используемые для получения точечного разложения Холецкого, для блочной матрицы [math]А[/math] останутся практически без изменений. Вместо явного обращения диагональных блоков, эффективнее будет хранить их в факторизованном виде [math]D_{ii}=L_{ii}L^T_{ii}[/math], а вместо точечной операции деления использовать операции решения для треугольных систем. Общее количество арифметических операций при этом останется практически неизменным, но зато существенно возрастет локальность вычислений. Размер блока [math]M[/math] выбирают таким образом, чтобы все блоки, участвующие в операции исключения, помещались в кэш первого или второго уровня. В этом случае подкачки данных в память будут минимальными.

Аналогичный прием понадобится также и для эффективной реализации параллельной версии разложения Холецкого, что позволит минимизировать как общее количество обменов, так и количество пересылаемой между процессорами информации. Полезным побочным эффектом использования блочной версии разложения Холецкого может стать повышение скалярной эффективности алгоритма за счет явного использования размера блока [math]M[/math] во внутренних циклах (прием "разворачивание цикла" или "loop unrolling").

3 Разложение Холецкого, точечный вещественный вариант для разреженной симметричной положительно-определённой матрицы

Если исходная матрица [math]A[/math] представлена в разреженном виде, то для экономии памяти, а также арифметических операций, необходимо учитывать ее разреженность.

3.1 Основные отличия от случая плотной матрицы

В этом разделе необходимо рассмотреть следующие виды разреженности:

1. ленточная матрица - матрица,ненулевые элементы которой сосредоточены внутри ленты шириной [math]2d+1[/math], т.е. когда [math]a_{ij}=0[/math] при [math]|i-j|\gt d[/math]. В этом случае, при проведении разложения Холецкого новые ненулевые элементы (если внутри ленты имеются нулевые элементы) могут образовываться только внутри этой же ленты. Количество ненулевых элементов в исходной матрице [math]A[/math], а также в нижнетреугольном множителе [math]L[/math] будет около [math](d+1)n[/math], а арифметические затраты соcтавят приблизительно [math]d^2n[/math]. Случай, когда матрица состоит всего из нескольких диагоналей (например, при конечно-разностной аппроксимации уравнений в частных производных на регулярной сетке) здесь отдельно не рассматривается, т.к. заполнение в нижнетреугольном множителе [math]L[/math] все-равно будет определяться исключительно наличием самой внешней (дальней) диагонали.

2. профильная матрица - в более общем случае, заполнение в каждой строке треугольного множителе [math]L[/math] будет определяться позицией первого ненулевого элемента. Сумма по всем строкам растояний от первого ненулевого элемента строки до главной диагонали и сотавляет "профиль" матрицы и определяет количество ненулевых элементов в нижнетреугольном множителе [math]L[/math].

3. матрица общей структуры разреженности. Верхней границей заполнения треугольного множителя [math]L[/math], конечно же, будет значение "профиля" матрицы, но учет особенностей структуры ненулевых элементов внутри профиля иногда может дать дополнительный эффект в повышении эффективности вычислений.

При рассмотрении общего случая разреженности необходимо выбрать формат хранения разреженных данных. Таковым может быть, например, формат построчного сжатия данных ("compressed sparse row" или CSR формат). В первом вещественном массиве, подряд (обычно в порядке возрастания номером столбцов) хранятся ненулевые элементы матрицы, во втором, в том же порядке хранятся номера столбцов, в третьем, отдельно сохраняется начало каждой строки. Если общее количество ненулевых элементов в матрице равно nnz ("number of nonzeros"), то память для хранения разреженных данных такой матрицы в формате CSR при использовании двойной точности составит [math]3\,{\rm nnz}+n+1[/math]. Оценку количества арифметических операций в общем случае невозможно, т.к. помимо количества ненулевых элементов в исходной матрице оно существенно зависит от структуры ее разреженности.

Для реализации разложения Холецкого в этом случае понадобится несколько операций с разреженными строками:

• копирование из одной разреженной строки в другую (или во временный "плотный" вектор, операция распаковки данных);

• выполнение операции исключения для одного из элементов строки;

• вставка в строку нового ненулевого элемента ("fill-in");

• сжатие данных с копированием из временного плотного вектора в сжатый разреженный (операция упаковки данных).

3.2 Переупорядочивания для уменьшения количества новых ненулевых элементов

Структура треугольного множителя [math]L[/math], а также объем памяти им занимаемый, зависят от упорядочивания строк и столбцов исходной матрицы [math]A[/math], в котором проводится разложение. Существуют алгоритмы, минимизирующие заполнение матрицы [math]L[/math].

• В первую очередь это алгоритм RCM (reversed Cuthill–McKee), который предназначен для уменьшения профиля матрицы. Одновременно с уменьшением профиля происходит и уменьшение заполнения треугольного множителя [math]L[/math]. Это очень широко применяемый, быстрый, но не самый эффективный алгоритм. Русского аналога название этого алгоритма не имеет.

• Алгоритм вложенных сечений (Nested Dissection, ND) - служит именно для минимизации заполнения множителя [math]L[/math]. В некоторых частных случаях доказана его ассимптотическая оптимальность.

В общем случае, проблема поиска перестановки минимизирующей заполнение множителя [math]L[/math] является NP-полной задачей.

4 Разложение Холецкого, блочный вещественный вариант для разреженной симметричной положительно-определённой матрицы

Иногда разреженную симметричную матрицу бывает удобно представить в блочном виде с блоками небольшого размера [math]M[/math], равного, например, количеству неизвестных функций на узел при конечно-элементной или конечно-разностной аппроксимации уравнений в частных производных. В этом случае структура разреженности хранится сразу для блочной структуры разреженности (что позволяет экономить память на хранении целочисленных массивов). Если общее количество ненулевых блоков размера [math]M\times M[/math] в матрице равно nnz ("number of nonzeros"), то пямять для хранения разреженных данных такой мелкоблочной матрицы в формате CSR при использовании двойной точности составит [math](2M^2+1)\,{\rm nnz}+n/M+1[/math].

В некоторых случаях, размер блока [math]M[/math] может выбираться искуственно, например, для повышения эффективности работы процедур нижнего уровня за счет приема разворачивания циклов (loop unrolling).

Алгоритмы, необходимые при выполнении разложения Холецкого для матриц, рассмотренных в этом разделе, могут быть получены комбинацией уже рассмотренных идей блочности и разреженности.

5 Разложение Холецкого для эрмитовой матрицы

Эрмитовой (или комплексно-самосопряженной) матрицей называют такую квадратную комплексную матрицу [math]A[/math], для элементов которой выполняется соотношение [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].

5.1 Точечный вариант

Как естественное обобщение разложения Холецкого для точечной симметричной положительно-определеной матрицы может быть рассмотрено разложение Холецкого для эрмитовой положительно-определеной матрицы. Все формулы для вычисления разложения остаются прежними, только теперь вместо операций над вещественными числами выполняются аналогичные комплексные операции:

[math] \begin{align} L_{ii} & = \sqrt{ A_{ii} - \sum_{k=1}^{i-1} L_{ik}L_{ik}^* }, \\ L_{ij} & = \frac{1}{L_{jj}} \left( A_{ij} - \sum_{k=1}^{j-1} L_{ik} L_{jk}^* \right), \quad j \lt i. \end{align} [/math]

В отличие от вещественного варианта, для выполнении аналогичных комплексных операций потребуется считывать из памяти вдвое больше данных и производить над ними примерно вчетверо больше арифметических операций, что должно не только несколько улучшить локальность вычислений, но и повысить их эффективность.

5.2 Блочный вариант

Реализация блочного варианта разложения Холецкого для эрмитовых матриц будет аналогична рассмотрему выше блочному варианту для вещественных матриц.

6 Использование разложения Холецкого в итерационных методах

При выполнении разложения Холецкого в арифметике с фиксированной машинной точностью полученные треугольный фактор [math]L[/math] и само решение может оказаться недостаточно точным. Для получения более точного решения может применяться некоторый итерационный метод (например, метод сопряженных градиентов), с использованием полученного разложения [math]LL^T[/math] в качестве предобуславливателя.

Основной причиной формирование неполного или неточного разложения в качестве предобуславливателя чаще всего бывает требование экономии памяти.

6.1 Ограничивание заполнения в разложении Холецкого

При выполнении разложения Холецкого для разреженной матрицы, может образовываться такое большое количество новых ненулевых элементов, что оперативной памяти на хранение полного разложения будет недостаточно. В этом случае можно построить неполное или приближенное разложение для применения его в дальнейшем в качестве предобуславливателя. В англоязычной литературе для обозначения таких разложения применяют единый термин Incomplete Cholesky factorization, или сокращенно IC разложение.

6.2 Неполное разложение Холецкого по позициям IC([math]k[/math])

Неполное разложение Холецкого можно получить используя заранеее выбранные ограничения по структуре заполнения. Чаще всего получают разложение Холецкого на тех же позициях, в которых находятся ненулевые элементы исходной матрицы [math]A[/math]. Такое разложение обозначают IC(0) или просто IC0.

Если качества разложения IC0 оказывается недостаточно, то можно выбрать более широкую структуру тругольного множителя [math]L[/math], например, разрешить образование одного уровня новых ненулевых элементов от исходной структуры матрицы [math]A[/math]. Формально, это означает заполнение внутри структуры матрицы [math]A^2[/math], а такое разложение обозначают IC(1).

Можно рассмотреть и более общий случай, с заполнением внутри структуры матрицы [math]A^{k+1}[/math], где [math]k \geq 0[/math]. Такое разложение обозначают IC([math]k[/math]).

Обычно с ростом значения [math]k[/math] точность неполного разложения IC([math]k[/math]) возрастает, хотя это совсем не является обязательным даже для симметричных положительно определенных матриц, полное разложение для которых существует и находится однозначно. Из-за неполноты разложения на главной диагонали могут оказаться нулевые или даже отрицательные элементы. Чтобы завершить треугольное разложение в таких случаях применяют предварительный диагональный сдвиг исходной матрицы [math]A+\varepsilon I[/math] перед ее разложением. Здесь [math]\varepsilon\gt 0[/math] - малый параметр, а [math]I[/math] - диагональная матрица. Если слишком малый или неположительный диагональный элемент образуется в процессе разложения, то применяют его замену на некоторое заранее выбранное значение. Такую операцию называют диагональной коррекцией разложения.

Неполное разложение IC([math]k[/math]) иногда называют также "разложение по позициям".

6.3 Приближенное разложение Холецкого по значениям IC([math]\tau[/math])

Для контроля заполнения в треугольном множителе [math]L[/math] разложения Холецкого, кроме структурных ограничений, можно также применить ограничение разложения в зависимости от значения самих элементов разложения. Например, можно сохранять только элементы, большие по модулю чем некоторый малый параметр [math]\tau\gt 0[/math]. В этом случае разложение называют приближенным разложением Холецкого или разложением "по значению" и обозначают IC([math]\tau[/math]). Величину [math]\tau[/math] называют "порогом" разложения или "порогом" фильтрации.

Вполне правомерным является ожидание того, что в уменьшением [math]\tau[/math] точность полученного разложения будет возрастать, правда за счет роста количества ненулевых элементов в треугольном множителе [math]L[/math]. Недостатком же такого разложения является то, что, в общем случае, предсказать заполнение [math]L[/math] не возможно.

С точки зрения устойчивости разложения вариант приближенного разложения Холецкого является более предпочтительным, хотя применение предварительного диагонального сдвига, а также диагональной коррекции также допускается. Если же описанные приемы не помогаю получить разложения достаточной точности, то можно применить прием модификации диагонали Азиза-Дженингса, который при отбрасывании малого элемента разложения [math]\ell_{ij}[/math] состоит в добавлении модуля этого элемента к диагональным элементам разложения [math]\ell_{ii}[/math] и [math]\ell_{jj}[/math]. Это прием гарантирует существование приближенного разложения для любой симметричной положительно определенной матрицы [math]A[/math]. Наиболее эффективно этот прием модификации главной диагонали можно организовать при использовании Ктаут-версии разложения Холецкого.

6.4 Приближенное разложение Холецкого второго порядка IC([math]\tau_1,\tau_2[/math])

Для повышения точности приближенного разложения можно применить "двухпороговую" версию приближенного разложения Холецкого. Основная идея такого разложения, назывемого разложением Тисменецкого-Капорина, состоит в том чтобы вычисление разложения проводить в более высокой точности [math]\tau_2[/math], а сохранять в треугольном множителе только значения, которые по модулю не меньше [math]\tau_1[/math]. Обычно полагают [math]\tau_1=\tau[/math] и [math]\tau_2=\tau^2[/math], в этом случае разложение называют разложением "второго порядка", т.к. элементы матрицы ошибок оказываются по модулю меньше чем [math]\tau^2[/math]. Для обозначения симметричного порогового разложения второго порядка используют обозначение IC2, которое не следует путать со структурным разложением IC(2) (т.е. разложением IC([math]k[/math]), где [math]k=2[/math]).

Такое разложение обычно используется вместе с приемом Азиза-Дженингса для модификации диагональных элементов, получая вариант "безотказного" разложения для любой симметричной положительно определенной матрицы [math]A[/math]. Этот вариант разложения позволяет получать наиболее точные разложения (при одинаковом заполнении множителя [math]L[/math]), хотя для их вычисления приходится тратить больше времени на вычисление самого разложения.

6.5 Комбинация разложений Холецкого IC([math]k,\tau[/math]) и IC([math]\tau,m[/math])

Для экономии памяти при вычислении неполного или приближенного разложения Холецкого можно использовать следующие два варианта симметричных треугольных разложений.

Для контроля верхней границы заполнения треугольного множителя [math]L[/math] можно предложить использовать заполнение как и для разложения IC([math]k[/math]), при некотором выбранном значении [math]k[/math]. Для дальнейшей экономии памяти разложение в заданной структуре разреженности можно вести с использованием порога разложения [math]\tau[/math], как и при проведении разложения IC([math]\tau[/math]). Такую комбинацию можно назвать IC([math]k,\tau[/math]) разложением. Применяться она может, например, при необходимости априорных структурных ограничений для минимизации обменов при использовании параллельной версии разложения для распределенной памяти.

Второй вариант структурно-порогового разложения можно описать следующим образом. При проведении обычного порогового IC([math]\tau[/math]) разложения, наложим дополнительное ограничение на элементы строк матрицы [math]L[/math]: разрешим сохранение только не более чем [math]m[/math] наибольших по модулю элементов рассматриваемой строки разложения [math]L[/math]. При общей размерности задачи [math]n[/math] в матрице [math]L[/math] будет не более чем [math]nm[/math] элементов. Такой подход представляется разумным, например, для матриц полученных в результате дискретизации с достаточно регулярным шаблоном. Наиболее известен несимметричный вариант такого разложения, предложенного Саадом и называемого разложением ILUT.

7 Использование разложения Холецкого в параллельных итерационных алгоритмах

Формулы разложения Холецкого по большей части имеют рекурсивный характер и выделение параллельных и независимых этапов вычислений является неочевидной и непростой задачей. Слишком прямолинейное ее решение может привести к значительному объему пересылаемых данных, что значительно снизит результат распараллеливания. Наибольший эффект может дать подход, основанный на предварительном переупорядочивании исходной матрицы.

7.1 Переупорядочивания для выделения блочности

Для того чтобы выделить независимые блоки вычислений, можно использовать симметричные перестановки строк и столбцов исходной матрицы, приводящие ее к блочно окаймленному виду. В этом случае основная часть работы будет сосредоточена в независимых блоках, которые можут обрабатываться параллельно и без обменов данными.

Наиболее простым, но не слишком эффективных способом упорядочивания является предварительное упорядочивание матрицы с помощью обратного алгоритма Катхилла-Макки (RCM, reverse Cuthill—McKee) для минимизации ширины профиля, а затем равномерное разбиение по блокам (процессорам) в новом упорядочивании. После присваивания номера процессора каждой вершине графа матрицы, в независимые блоки можно выделить те вершины графа, которые связаны только с вершинами, имеющими тот же номер процессора (т.е. являющимися внутренними вершинами). Остальные вершины можно объединить в последний блок окаймления, который будет обрабатываться отдельно. Все вычисления внутри блоков будут независимы и могут выполняться параллельно. Для повышения эффективности треугольной факторизации внутренние вершины каждого из блоком можно также упорядочить с помощью метода RCM.

Более эффективными с точки зрения минимизации ширины окаймления будут следущие методы:

• Метод минимальных сепараторов, который заключается в последовательном нахождении имеющих минимальный размер разделителей (сепараторов), обеспечивающих расщепление оставшихся вершин на два независимых блока.

• Метод минимальной степени (Minimum Degree, MD). Прямое применение этого метода к матрицах большого размера затруднительно, поэтому используется приближенный метод минимальной степени (Approximate Minimum Degree, AMD).

• Метод вложенных сечений (Nested Dissection, ND). Это рекурсивный алгоритм, на каждом шаге разделяющий множество вершин на два нерависимых блока, представленые в [math]2\times2[/math] блочно-окаймленном виде.

В качестве побочного положительного эффекта от такого упорядочивания, для некоторого вида матриц доказано, что полное разложение будет иметь почти минимально возное количество ненулевых элементов. Нахождение оптимальной перестановки в общем случае является NP-полной задачей.

Существуют и другие алгоритмы упорядочивания матриц для наиболее оптимального их распределения по процессорам. Наиболее популярными являются последовательные пакеты METIS, JOSTLE, SCOTCH, CHACO, PARTY, а также параллельные коды PARMETIS, JOSTLE, PT-SCOTCH и ZOLTAN. Многие из них являются свободно распространяемыми.

7.2 Разложение в независимых блоках

Вычисления в независимых блоках полностью независимы и могут выполняться паралельно без обменов данными между процессорами. Единственным недостатком может быть лишь то, что для пороговых разложений количество арифметических опреаций на различных процессорах может различаться, что может привести к некоторой несбалансированности вычислений.

7.3 Разложение в сепараторах

Последний блок, в котором собраны сепараторы всех блоков при небольшом количестве используемых процессоров может обрабатываться, например, на одном головном процессоре. Если же процессоров достаточно много, но обработку сепараторов также необходимо производить совместно.

7.4 Иерархические и вложенные алгоритмы

Для обработки сепараторов могут быть применены те же алгоритмы упорядочивания и выделения независимых блоков что и для исходной задачи. Этот же прием может быть применен и на следующем шаге, с получением иерархического или вложенного алгоритма.

В случае применения порогового разложения такие построения могут быть явно применены к построенному дополнению по Шуру. Для структурных факторизаций существуют приложения сразу строящие многоуровневые упорядочивания и обеспечивающие минимальность заполнения и обменов.

7.5 Блочный метод Якоби

Описанные выше подходы относились к типу "прямого" или "явного" разложения, когда при отбрасывании элементов разложения в расчет принимались только их абсолютные значения. Структурные свойства при этом являлись подчиненными и на отбрасывание элементов никак не влияли.

Алтернативням является подход, когда из-за соображений увеличения ресурса параллелизма некоторые элементы отбрасываются исключительно из структурных соображений. Например, можно сначала распределить строки матрицы по процессорам, а затем перед проведением разложения просто отбросить все элементы связывающие один процессор с другим. В этом случае разложение будет проходить полностью независимо в каждом из блоков. Внутри каждого блока может проводиться любой из структурных или пороговых вариантов разложения Холецкого. Построение такого предобуславливателя называют блочным методом Якоби без перекрытия блоков (Block Jacobi, BJ). Такое предославливание является наиболее простым, применение его наиболее параллельным (полностью без обменов данными), правда сходимость может оставлять желать лучшего, почти независимо от качества разложения внутри каждого из блоков.

7.6 Адитивный метод Шварца

разложение гораздо более высокого качества по сравнению с мотодом Якоби можно получить применяя адитивный метод Шварца (Additive Schwarz, AS). Иногда этот метод называют также методом Якоби с перекрытиями. Суть его заключается в расширении структуры каждого из блоков матрицы за счет добавления нескольких слоев близлежащих строк матрицы. Треугольное разложение строится для расширенной матрицы, таким образом на каждом из процессоров происходит решение расширенной подзадачи с привлечением данных от других процессоров. После нахождения решения подзадачи на каждом из процессоров обычно происходит отбрасывание не локальных компонент решения. Такой вариант метода называют адитивный метод Шварца с ограничениями (Restricted Additive Schwarz, RAS).

Сходимость адитивного метода Шварца бывает гораздо выше чем сходимость метода Якоби, и обычно монотонно улучшается с ростом размера перекрытия. Несмотря на дополнительные вычисление и обмены, общее время решения на параллельном компьютере может быть существенно меньше.

7.7 Неполное обратное треугольное разложения

Существует и другой вариант аддитивного разложения, который кроме несколько более быстрой сходимости опирается на построение перекрытий блоков только в одну сторону ("назад", т.е. в сторону меньших номеров строк). Название этого метода блочное неполного обратного разложения Холецкого, имеющее только английскую аббревиатуру BIIC (Block Incomplete Inverse Cholesky). Позднее, вместе с рассмотрением несимметричного варианта разложения, этот метод стал называться методом неполного обратного треугольного разложения, или НОТ разложения.

Комбинация этого метода с неполным симметричным треугольным разложением второго порядка IC2 внутри каждого из блоков имеет обозначение BIIC2.

Идея этого метода впервые была предложена И.Е.Капориным без рассмотрения возможностей его распараллеливания. В литературе встречается также название этого метода как метод Капорина-Коньшина, по имени авторов впервые представивших его параллельную реализацию и проанализировавших ее свойства.

8 Решение линейных систем с треугольной матрицей

Разложение Холецкого может применяться для решения системы линейных уравнений [math]Ax = b[/math], если матрица [math]A[/math] симметрична и положительно-определена. Выполнив разложение [math]A = LL^T[/math], решение [math]x[/math] получается последовательным решением двух треугольных систем уравнений [math]Ly = b[/math] и [math]L^T x = y[/math].

8.1 Решение системы с плотной нижнетреугольной матрицей

Решение линейной системы с плотной нижнетреугольной матрицей [math]L y = b[/math] можно представить в виде "прямого" хода, т.е. цепочки вычислений, начиная с верхнего угла матрицы [math]L[/math] по возрастанию номера строки [math]i[/math]:

[math] \begin{align} y_{1} & = b_{1}, \\ y_{i} & = b_{i} - \sum_{j = 1}^{i-1} \ell_{ij} y_{j}, \quad i = 2,...,n. \end{align} [/math]

В разделе Прямая_подстановка_(вещественный_вариант) содержится подробное описание алгоритма и его анализ.

8.2 Решение системы с плотной верхнетреугольной матрицей

Решение линейной системы с плотной верхнетреугольной матрицей [math]U x = y[/math] (где, например, [math]U=L^T[/math]) можно представить в виде "обратного" хода, т.е. цепочки вычислений, начиная с нижнего угла матрицы [math]U[/math] при убываниии номера строки [math]i[/math]:

[math] \begin{align} x_{n} & = y_{n}/u_{nn}, \\ x_{i} & = \left (y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} \right ) / u_{ii}, \quad i = n - 1,...,1. \end{align} [/math]

В разделе Обратная_подстановка_(вещественный_вариант) содержится подробное описание алгоритма и его анализ.

8.3 Решение системы с разреженной нижнетреугольной матрицей

Решение линейных систем с разреженной нижне- или верхнетреугольной матрицей аналогично рассмотренным алгоритмам для плотных матриц, при этом подстановки ведутся исключительно для ненулевых элементов с учетом идеи работы с разреженными матрицами.

8.4 Решение системы с комплексной треугольной матрицей

Решение линейных систем с комплексной нижне- или верхнетреугольной матрицей аналогично рассмотренным алгоритмам для вещественных матриц, при этом арифметические операции выполняются в комплексной арифметике, аналогично операциям раздела факторизации эрмитовых матриц.

8.5 Решение систем с блочноокаймленными треугольными матрицами

Особенность решения линейных систем с блочноокаймленными треугольными матрицами в том что независимость вычислений в отдельных блоках дает возможность проведения параллельных вычислений.

9 Существующие реализации алгоритма

• В LAPACK используется функция DPBTRF (последовательная реализация для двойной точности).

• В ScaLAPACK используется функция PDPBTRF (паралельная реализация для двойной точности).

• В SAS используется функция ROOT( matrix ), входящая в пакет SAS IML.

• В системах MATLAB, Octave, R разложение выполняется командой U = chol(A).

• В Maple и NumPy существует процедура cholesky в модуле linalg.

• В Mathematica используется процедура CholeskyDecomposition[A].

• В GSL используется функция gsl_linalg_cholesky_decomp.

• В ALGLIB имеются реализации как LLT так и LDLT разложений для различных языков програмирования: C#, C++, C++ (арифметика повышенной точности), FreePascal, Delphi, VB.NET, VBA, Python.

• В Online Matrix Calculator непосредственно в web-интерфейсе можно выполнить разложение Холецкого, выбрав раздел Cholesky Decomposition.