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

Материал из Алговики
Перейти к навигации Перейти к поиску
[непроверенная версия][непроверенная версия]
м (/* Разложение Холецкого (метод квадратного корня), базовый точечный вещественный вариант для плотной симметричной положительно-опреде…)
Строка 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>A</math> в виде произведения <math>A = LL^T</math>, где <math>L</math> — нижняя (''Lower'') треугольная матрица со строго положительными элементами на диагонали. Иногда разложение удобно записать в эквивалентной форме <math>A = U^TU</math>, где <math>U = L^T</math> — верхняя (''Upper'') треугольная матрица. Разложение Холецкого всегда существует и единственно для любой симметричной положительно-определённой матрицы.
Строка 21: Строка 23:
 
В разделе [[Разложение Холецкого (метод квадратного корня)]] подробно рассмотрен базовый точечный вещественный вариант для плотной симметричной положительно-определённой матрицы.
 
В разделе [[Разложение Холецкого (метод квадратного корня)]] подробно рассмотрен базовый точечный вещественный вариант для плотной симметричной положительно-определённой матрицы.
  
Иногда удобнее бывает рассматривать <math>LDL^T</math> вариант разложения Холецкого, в котором матрица <math>L</math> имеет единицы на главной диагонали, а <math>D</math> - диагональная матрица с положительными элементами. В этом варианте разложения легко проследить связь с ранее рассмотренным <math>LL^T</math> вариантом:
+
=== <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>
 
<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>

Версия 13:19, 19 марта 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] \ell_{ii} = \sqrt{a_{ii} - \sum_{k=1}^{i-1} \ell_{ik}^2} [/math]

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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 Решение системы с плотной верхнетреугольной матрицей

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

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

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

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

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

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.