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

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

[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{ \ell_{ii} = \sqrt{a_{ii} - \sum_{k=1}^{i-1} \ell_{ik}^2} }[/math]

[math]\displaystyle{ \ell_{ij} = \frac{1}{\ell_{jj}} \left(a_{ij} - \sum_{k=1}^{j-1} \ell_{ik} \ell_{jk} \right), ~если~ j \lt i. }[/math]

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

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

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

2) Краут-вариант алгоритма Холецкого (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{ L_{ii} = \sqrt{ A_{ii} - \sum_{k=1}^{i-1} L_{ik}L_{ik}^* } }[/math]

[math]\displaystyle{ L_{ij} = \frac{1}{L_{jj}} \left( A_{ij} - \sum_{k=1}^{j-1} L_{ik} L_{jk}^* \right), ~если~ j \lt i. }[/math]

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

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

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

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

Разложение Холецкого может применяться для решения системы линейных уравнений [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].

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

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

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