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

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

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

1.1 $LL^T$-разложение

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

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

$\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}$

Выражение под квадратным корнем всегда положительно, если $A$ — вещественная симметричная положительно определённая матрица.

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

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

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

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

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

1.2 $LDL^T$-разложение

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

$A = LDL^T = LD^{1/2}D^{1/2}L^T = (LD^{1/2})\,(LD^{1/2})^T = \tilde L \tilde L^T,$

так и с несимметричным $LU$-разложением:

$A = LDL^T = L(DL^T) = LU.$

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

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

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

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

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

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

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

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

• Профиль - в более общем случае, заполнение в каждой строке треугольного множителе $L$ будет определяться позицией первого ненулевого элемента. Сумма по всем строкам расстояний от первого ненулевого элемента строки до главной диагонали и составляет "профиль" матрицы и определяет верхнюю границу количества ненулевых элементов в нижнетреугольном множителе $L$.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 Разложение Холецкого для симметричной незнакоопределенной (седловой) матрицы

Если симметричная матрица ${\mathcal A}$ представима в виде

${\mathcal A} = \begin{bmatrix} A & B^T \\ B & -C \end{bmatrix} ,$

где $A$ - симметричная положительно определенная ($A=A^T\gt 0$) и $C$ - симметричная неотрицательно определенная ($C=C^T\ge0$) матрицы, то, выполнив один шаг блочного исключения, ее можно преобразовать к виду

$\begin{bmatrix} A & 0 \\ 0 & S \end{bmatrix} ,$

где матрица дополнения по Шуру $S=-(C+B^TA^{-1}B)$ является строго отрицательно определенной ($S=S^T\lt 0$). Это означает, что матрица ${\mathcal A}$ имеет $n_A$ положительных и $n_C$ отрицательных собственных значений, где через $n_A$ и $n_C$ обозначены размерности матриц $A$ и $C$, соответственно.

В этом случае существует симметричное треугольное разложение вида ${\mathcal A}={\mathcal L}D{\mathcal L}^T$, где ${\mathcal L}$ является нижней унитреугольной, а диагональная матрица $D$ содержит $n_A$ положительных и $n_C$ отрицательных элементов на главной диагонали, причем такое разложение может быть получено напрямую без выбора ведущего элемента даже если $C$ - нулевая матрица.

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

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

Эрмитовой (или комплексно-самосопряженной) матрицей называют такую квадратную комплексную матрицу $A$, для элементов которой выполняется соотношение $a_{ij}=\overline{a_{ji}}$ (здесь, если $z=a+{\rm i\,}b\,$ и ${\rm i}^2=-1$, то $\overline z=a-{\rm i\,}b\,$). В матричном виде это можно записать как $A=\overline{A^T}$ или $A=A^*=A^Н$.

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

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

$\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}$

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

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

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

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

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

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

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

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

7.2 Неполное разложение Холецкого по позициям IC($k$)

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

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

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

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

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

7.3 Приближенное разложение Холецкого по значениям IC($\tau$)

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

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

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

7.4 Приближенное разложение Холецкого второго порядка IC($\tau_1,\tau_2$)

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

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

7.5 Комбинация разложений Холецкого IC($k,\tau$) и IC($\tau,m$)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

$\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}$

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

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

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

$\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}$

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

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

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

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

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

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

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

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

• В 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.