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

Материал из Алговики
Перейти к навигации Перейти к поиску
[непроверенная версия][непроверенная версия]
Строка 119: Строка 119:
  
 
== Использование разложения Холецкого в итерационных методах ==
 
== Использование разложения Холецкого в итерационных методах ==
 +
 +
При выполнении разложения Холецкого в арифметике с фиксированной машинной точностью полученные треугольный фактор <math>L</math> и само решение может оказаться недостаточно точным.
 +
Для получения более точного решения может применяться некоторый итерационный метод (например, метод сопряженных градиентов), с использованием полученного разложения <math>LL^T</math> в качестве предобуславливателя.
 +
 +
Основной причиной формирование неполного или неточного разложения в качестве предобуславливателя чаще всего бывает требование экономии памяти.
  
 
=== Ограничивание заполнения в разложении Холецкого ===
 
=== Ограничивание заполнения в разложении Холецкого ===
 +
 +
При выполнении разложения Холецкого для разреженной матрицы, может образовываться такое большое количество новых ненулевых элементов, что оперативной памяти на хранение полного разложения будет недостаточно.
 +
В этом случае можно построить неполное или приближенное разложение для применения его в дальнейшем в качестве предобуславливателя.
 +
В англоязычной литературе для обозначения таких разложения применяют единый термин Incomplete Cholesky factorization, или сокращенно IC разложение.
  
 
=== Неполное разложение Холецкого по позициям IC(<math>k</math>) ===
 
=== Неполное разложение Холецкого по позициям 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>0</math> - малый параметр, а <math>I</math> - диагональная матрица.
 +
Если слишком малый или неположительный диагональный элемент образуется в процессе разложения, то применяют его замену на некоторое заранее выбранное значение. Такую операцию называют диагональной коррекцией разложения.
 +
 +
Неполное разложение IC(<math>k</math>) иногда называют также "разложение по позициям".
  
 
=== Приближенное разложение Холецкого по значениям IC(<math>\tau</math>) ===
 
=== Приближенное разложение Холецкого по значениям IC(<math>\tau</math>) ===
 +
 +
Для контроля заполнения в треугольном множителе <math>L</math> разложения Холецкого, кроме структурных ограничений, можно также применить ограничение разложения в зависимости от значения самих элементов разложения.
 +
Например, можно сохранять только элементы, большие по модулю чем некоторый малый параметр <math>\tau>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>. Наиболее эффективно этот прием модификации главной диагонали можно организовать при использовании Ктаут-версии разложения Холецкого.
  
 
=== Приближенное разложение Холецкого второго порядка IC(<math>\tau_1,\tau_2</math>) ===
 
=== Приближенное разложение Холецкого второго порядка 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>.
 +
 +
Такое разложение обычно применяют вместе с приемом Азиза-Дженингса для модификации диагональных элементов, получая вариант "безотказного" разложения для любой симметричной положительно определенной матрицы <math>A</math>.
 +
Этот вариант разложения получает получать наиболее точные разложения (при одинаковом заполнении множителя <math>L</math>), хотя для их вычисления приходится тратить больше времени на вычисление самого разложения.
  
 
=== Комбинация разложений Холецкого IC(<math>k,\tau</math>) и IC(<math>\tau,m</math>) ===
 
=== Комбинация разложений Холецкого 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.
  
 
== Использование разложения Холецкого в параллельных итерационных алгоритмах ==
 
== Использование разложения Холецкого в параллельных итерационных алгоритмах ==

Версия 15:48, 25 марта 2015

Содержание

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]

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

Такое разложение обычно применяют вместе с приемом Азиза-Дженингса для модификации диагональных элементов, получая вариант "безотказного" разложения для любой симметричной положительно определенной матрицы [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 Переупорядочивания для выделения блочности

7.1.1 Метод минимальных сепараторов

7.1.2 Метод минимальной степени (Minimum Degree - MD)

7.1.3 Метод вложенных сечений (Nested Dissection - ND)

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

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

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

7.5 Блочный метод Якоби (без перекрытия блоков, Block Jacobi - BJ)

7.6 Адитивный метод Шварца (Additive Schwarz - AS)

7.7 Блочный метод неполного обратного разложения Холецкого (BIIC)

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.