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

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.

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.