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

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

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

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

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

$\begin{align} \ell_{ii} & = \sqrt{a_{ii} - \sum_{k=1}^{i-1} \ell_{ik}^2}, \\ \ell_{ij} & = \frac{1}{\ell_{jj}} \left(a_{ij} - \sum_{k=1}^{j-1} \ell_{ik} \ell_{jk} \right), \quad j \lt i. \end{align}$

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

$\begin{align} L_{ii} & = \sqrt{ A_{ii} - \sum_{k=1}^{i-1} L_{ik}L_{ik}^* }, \\ L_{ij} & = \frac{1}{L_{jj}} \left( A_{ij} - \sum_{k=1}^{j-1} L_{ik} L_{jk}^* \right), \quad j \lt i. \end{align}$

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

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

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

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

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

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

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

$\begin{align} y_{1} & = b_{1}, \\ y_{i} & = b_{i} - \sum_{j = 1}^{i-1} \ell_{ij} y_{j}, \quad i = 2,...,n. \end{align}$

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

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

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

$\begin{align} x_{n} & = y_{n}/u_{nn}, \\ x_{i} & = \left (y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} \right ) / u_{ii}, \quad i = n - 1,...,1. \end{align}$

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

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.

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