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

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

[math]\displaystyle{ LL^T }[/math] разложение

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

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

[math]\displaystyle{ \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]\displaystyle{ A }[/math] — действительная положительно-определённая матрица.

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

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

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

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

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

[math]\displaystyle{ LDL^T }[/math] разложение

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[math]\displaystyle{ \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]

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[math]\displaystyle{ \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]

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

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

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

[math]\displaystyle{ \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]

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

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

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

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

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

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

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

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

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