Cholesky method

1 The Cholesky decomposition (or the square-root method): a basic dot version for dense real symmetric positive definite matrices

1.1 The [math]LL^T[/math] decomposition

The Cholesky decomposition (or the Cholesky factorization) is a decomposition of a symmetric positive definite matrix [math]A[/math] into the product [math]A = LL^T[/math], where the factor [math]L[/math] is a lower triangular matrix with strictly positive diagonal elements. In some cases it is convenient to rewrite this decomposition in its equivalent form [math]A = U^TU[/math], where [math]U = L^T[/math] is an upper triangular matrix. The Cholesky decomposition always exists and is unique for any symmetric positive definite matrix, since the diagonal elements of the matrix [math]L[/math] are strictly positive.

The following formulas should be used to find the elements of the matrix [math]L[/math] starting with the upper left-hand corner of the original matrix [math]A[/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]

The above expression under the square root sign is always positive if the matrix [math]A[/math] is real symmetric and positive definite.

The computations are performed from the top down and from left to right: the value of [math]L_{ij}[/math] ([math]j \lt i[/math]) is first determined and, then, the value of [math]L_{ii}[/math] is found. The following computational algorithms are usually employed:

• the Cholesky–Banachiewicz algorithm (or the Cholesky algorithm) when the computations are started from the upper left-hand corner of the matrix [math]L[/math] in the row-wise manner; this version of the decomposition is most commonly used, especially when the row-wise format of matrix storage is chosen for the matrix [math]L[/math];

• the Crout version of the Cholesky algorithm (or Cholesky–Crout algorithm) when the computations are also started from the upper left-hand corner of the matrix [math]L[/math] but in a column-wise order; this version is more rarely used and is employed when the column-wise format of matrix storage is chosen for the matrix [math]L[/math] and when the pivot elements should be modified in approximate decompositions.

Both these algorithms can be used when the factor [math]L[/math] is stored in the lower triangular part of the original matrix [math]A[/math].

A basic dot version of the Cholesky algorithm for dense real symmetric positive definite matrices is extensively analyzed in The Cholesky decomposition (the square root method).

1.2 The [math]LDL^T[/math] decomposition

In a number of cases, it is more convenient to use the [math]LDL^T[/math] version of the symmetric triangular decomposition when the matrix [math]L[/math] is lower unitriangular (its diagonal elements are equal to 1) and [math]D[/math] is a diagonal matrix with positive elements. This version of decomposition can easily be related to the above [math]LL^T[/math] version:

[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]

and with the [math]LU[/math] decomposition:

[math]A = LDL^T = L(DL^T) = LU[/math].

2 The Cholesky decomposition: a block version for dense real symmetric positive definite matrices

A block version of the Cholesky decomposition can also be considered. Let [math]n=MN[/math]. Then, the original matrix [math]A[/math] of order [math]n[/math] can be represented as a block matrix of order [math]N[/math] with blocks of order [math]M\lt math\gt . All the formulas used in the dot version of the Cholesky decomposition remain almost the same for the block representation of the matrix \lt math\gt А[/math]. Instead of inverting the diagonal blocks, it is more efficient to store these blocks in their factored form [math]D_{ii}=L_{ii}L^T_{ii}[/math] and to solve the corresponding linear triangular systems. The total number of arithmetic operations remains almost the same, but the locality of the algorithm increases considerably. The block order [math]M[/math] should be chosen in such a way that all the blocks used during elimination can be stored in the first-level or second-level cache. Hence, the data swapping becomes minimal.

A similar technique is also useful for an efficient implementation of a parallel version for the Cholesky decomposition, which allows one to minimize the total number of data exchange and the amount of data transferred among the processors. Another advantage of employing the block version of the Cholesky decomposition consists in the fact that the scalar efficiency of the algorithm can be increased by using the loop unrolling technique for the inner loops over the block order [math]M[/math].

3 The Cholesky decomposition: a dot version for real sparse symmetric positive definite matrices

If the original matrix [math]A[/math] is sparse, it is reasonable to take into account this property to reduce the storage an computational cost.

3.1 Major distinctions from the case of dense matrices

Here we consider the following sparsity types that characterize the ways of matrix storage.

1. Band. If the nonzero elements of a symmetric matrix are concentrated inside a band of width [math]2d+1[/math], then this matrix is called a symmetric band matrix. More precisely, a symmetric matrix [math]A[/math] is said to be a symmetric band matrix if [math]a_{ij}=0[/math] for [math]|i-j|\gt d[/math]. During the Cholesky decomposition, new nonzero elements that originally were zero may appear only inside this band, since no pivoting is done because of positive definiteness. In the original matrix [math]A[/math] and in the lower triangular factor [math]L[/math], the number of nonzero elements is about [math](d+1)n[/math], whereas the number of arithmetic operations is about [math]d^2n[/math]. Note that the situation when an original zero element becomes nonzero during elimination is also known as the generation of fill-in elements (or the generation of fill-ins).

2. Profile. In a more general case. the generation of fill-ins in each row of the lower triangular factor [math]L[/math] is specified by the position of the first nonzero element. The sum of the distances from the first elements of all rows to the main diagonal forms the "profile" of a matrix. In other words, the profile of a symmetric matrix is a measure of how close its elements are to the main diagonal and specifies an upper bound for the number of nonzero elements in the factor [math]L[/math].

3. General sparsity pattern. An upper bound for the number of nonzero elements in the factor [math]L[/math] is the value of its "profile"; however, the consideration of the peculiarities in the structure of nonzero elements inside the profile may sometimes increase the performance efficiency.

In the case of general sparse matrices, it is necessary to choose a storage format for matrix elements. As an example, we mention the "compressed sparse row" (CSR) format. In this format, the nonzero elements are stored row by row in the first real array; the second integer array contains the column indices of the nonzero elements stored in the first array; and the third integer array contains the pointers (or indices) to the beginning of each row in the first and second arrays. By nnz ("number of nonzeros") we denote the total number of the nonzero elements. In double precision, then, the amount of storage required by the CSR format is equal to [math]3\,{\rm nnz}+n+1\lt math\gt , where \lt math\gt n\lt math\gt is the number of rows in the matrix. Generally, the amount of arithmetic operations cannot be estimated, since this amount depends not only on the number of nonzero elements but is essentially dependent on the sparsity structure of the matrix. In order to implement the Cholesky decomposition in this case, the following operations should be used for sparse rows. * ''copying'' a sparse row to another one (or to a working "dense" vector, the data ''unpacking'' operation); * ''elimination operation'' for one of the row elements; * ''insertion'' of a new nonzero element (known also as a "fill-in" element) into a row; * ''data compression'' with copying a working dense vector to a compressed sparse vector (the data ''packing'' operation). === Reordering to reduce the number of fill-in elements === The structure of the lower triangular factor \lt math\gt L[/math] and the amount of the required storage are dependent on the ordering of the rows and columns of the original matrix [math]A[/math]. There exist algorithms minimizing the number of fill-ins in the matrix [math]L[/math]. Here we mention the following algorithms of such a type.

• The RCM (Reverse Cuthill–McKee) algorithm to reduce the profile of a sparse symmetric matrix (or its bandwidth). In addition to the profile reduction, it is possible to reduce the number of fill-ins in the factor [math]L[/math]. This algorithm is widely use in practice but is not the most efficient.

• The ND (Nested Dissection) algorithm is intended to minimize the number of fill-ins in the factor [math]L[/math]. Its asymptotic optimality can be proved for some particular cases.

In the general case, the problem of searching for the permutation that minimizes the number of fill-ins in the factor [math]L[/math] can be considered as an NP-complete problem.

4 The Cholesky decomposition: a block version for sparse real symmetric positive definite matrices

In a number of cases, it is convenient to represent a sparse symmetric matrix in a block form with blocks of small order [math]M[/math]. For example, the value of [math]M[/math] can be chosen to be equal to the number of unknown functions per a node when using various finite element or finite difference approximations of partial differential equations. The sparsity structure is stored in computer memory for the entire block sparsity structure, which allows one to save storage for integer arrays. Let nnz be the number of nonzero blocks of order [math]M[/math]. Then, the amount of memory storage required for such small block matrices is equal to [math](2M^2+1)\,{\rm nnz}+n/M+1[/math] if the CSR format and the double precision mode are used.

In some cases a value of the block order [math]M[/math] can be chosen from other reasonings (for example, on the basis of loop unrolling technique to increase the performance efficiency of low-level procedures).

Other algorithms devoted to the Cholesky decomposition for the matrices considered in the above discussion can be developed on the basis of the above approaches for block and [[Метод_Холецкого_(нахождение_симметричного_треугольного_разложения)#Разложение Холецкого, точечный вещественный вариант для разреженной симметричной положительно-определённой матрицы|sparse] matrices.

5 The Cholesky decomposition for symmetric indefinite (saddle point) matrices

If a symmetric matrix [math]{\mathcal A}[/math] is represented as

[math] {\mathcal A} = \begin{bmatrix} A & B^T \\ B & -C \end{bmatrix}, [/math]

where [math]A[/math] is a symmetric positive definite matrix ([math]A=A^T\gt 0[/math]) and [math]C[/math] is a symmetric nonnegative definite matrix ([math]C=C^T\ge0[/math]), then the original matrix can be transformed to the following form by performing a single step of block pivoting:

[math] \begin{bmatrix} A & 0 \\ 0 & S \end{bmatrix}. [/math]

Here the Schur complement matrix [math]S=-(C+B^TA^{-1}B)[/math] is a strictly negative definite matrix ([math]S=S^T\lt 0[/math]). This means that the matrix [math]{\mathcal A}[/math] has [math]n_A[/math] positive and [math]n_C[/math] negative eigenvalues, where [math]n_A[/math] and [math]n_C[/math] are the orders of the matrices [math]A[/math] and [math]C[/math], respectively.

In this case there exists a symmetric triangular decomposition [math]{\mathcal A}={\mathcal L}D{\mathcal L}^T[/math] such that [math]{\mathcal L}[/math] is a lower unitriangular matrix and the diagonal matrix [math]D[/math] consists of [math]n_A[/math] positive and [math]n_C[/math] negative elements. Such a decomposition can be obtained without pivoting, even if [math]C[/math] is a zero matrix.

In the general case of a symmetric nonsingular indefinite matrix, its decomposition should be performed by diagonal pivoting, which corresponds to a symmetric permutation of the rows and columns in the original matrix [math]{\mathcal A}[/math].

6 The Cholesky decomposition for Hermitian matrices

A Hermitian matrix is a square complex matrix [math]A[/math] such that the following equalities are valid for its elements: [math]a_{ij}=\overline{a_{ji}}[/math]; here we use the complex conjugation operation defined as follows: if [math]z=a+{\rm i\,}b\,[/math] and [math]{\rm i}^2=-1[/math], then [math]\overline z=a-{\rm i\,}b\,[/math]). In the matrix form, the above equalities can be written as [math]A=\overline{A^T}[/math] or [math]A=A^*=A^Н[/math]. Thus, a Hermitian matrix is equal to its conjugate transpose. Its diagonal elements are real.

6.1 A dot version

A dot version of the Cholesky decomposition for real symmetric positive definite matrices can be generalized to the case of Hermitian positive definite matrices. The above formulas derived for the real case remain the same with the only exception that the complex operations are used instead of the corresponding operations with real numbers:

[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]

Another distinction of this version from the real version consists in the following: for the complex operations, it is required to read a double amount of data from memory and to perform about four times as many arithmetic operations, which increases the locality of computations and the total performance efficiency.

6.2 A block version

An implementation of a block version of the Cholesky decomposition for Hermitian matrices is similar to the case of the block versions for real matrices.

7 Usage of Cholesky decomposition in iterative methods

The triangular factor [math]L[/math] and a solution to a corresponding linear system may not be accurate enough because of machine arithmetic. In order to refine the solution, a number of iterative methods (for example, the conjugate gradient method) can be employed using the [math]LL^T[/math] decomposition as a preconditioner.

The memory saving is the main reason to use an incomplete or inaccurate decomposition as a preconditioner.

7.1 The fill-in reduction in the Cholesky decomposition

When the Cholesky decomposition is used for sparse matrices, a large number of nonzero elements may appear. As a result, the complete decomposition may not be stored in the main memory. In this case it is reasonable to compute an incomplete or approximate decomposition for employing it as a preconditioner.

7.2 IC([math]k[/math]: the incomplete Cholesky decomposition (or factorization) by position

An incomplete Cholesky decomposition can be obtained using a chosen strategy of dropping fill-in elements while computing the incomplete factors. More often, the incomplete decomposition is stored in the same positions where the nonzero elements of the original matrix [math]A[/math] are stored. Such a decomposition is denoted by IC(0) or IC0.

If the quality of the IC0 decomposition is not enough, then it is reasonable to choose a wider structure of the triangular factor [math]L[/math]. For example, we can allow the formation of another level of new nonzero elements compared to the original structure of the matrix [math]A[/math]. Formally, in this example the generation of fill-in elements is allowed within the sparsity pattern of the matrix [math]A^2[/math]. Such a decomposition is denoted by IC(1).

It is also reasonable to consider a more general case when the generation of fill-in elements is allowed within the sparsity pattern of the matrix [math]A^{k+1}[/math], where [math]k \geq 0[/math]. Such a decomposition is denoted by IC([math]k[/math]).

Usually, when the value of [math]k[/math] increases, the accuracy of the IC([math]k[/math]) incomplete decomposition also increases; however, this is not always the case even for symmetric positive definite matrices, although their complete decompositions exist and are unique. Since the IC([math]k[/math]) decomposition is not complete, zero or negative elements may appear on the main diagonals. In these cases the diagonal shift [math]A+\varepsilon I[/math] is preliminary performed in an original matrix [math]A[/math] before its decomposition (here [math]\varepsilon\gt 0[/math] is a small parameter and [math]I[/math] is a diagonal matrix). If, during the decomposition, a too small or nonpositive element is formed on the main diagonal, then such an element is replaced by a prescribed value. This operation is called a diagonal correction of the decomposition.

7.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]. Наиболее эффективно этот прием модификации главной диагонали можно организовать при использовании Ктаут-версии разложения Холецкого.

7.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]. Для обозначения симметричного порогового разложения второго порядка используют обозначение IC2, которое не следует путать со структурным разложением IC(2) (т.е. разложением IC([math]k[/math]), где [math]k=2[/math]).

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

7.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 Использование разложения Холецкого в параллельных итерационных алгоритмах

Формулы разложения Холецкого по большей части имеют рекурсивный характер и выделение параллельных и независимых этапов вычислений является неочевидной и непростой задачей. Слишком прямолинейное ее решение может привести к значительному объему пересылаемых данных, что значительно снизит результат распараллеливания. Наибольший эффект может дать подход, основанный на предварительном переупорядочивании исходной матрицы.

8.1 Переупорядочивания для выделения блочности

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

Наиболее простым, но не слишком эффективных способом упорядочивания является предварительное упорядочивание матрицы с помощью обратного алгоритма Катхилла-Макки (RCM, reverse Cuthill—McKee) для минимизации ширины профиля, а затем равномерное разбиение по блокам (процессорам) в новом упорядочивании. После присваивания номера процессора каждой вершине графа матрицы, в независимые блоки можно выделить те вершины графа, которые связаны только с вершинами, имеющими тот же номер процессора (т.е. являющимися внутренними вершинами). Остальные вершины можно объединить в последний блок окаймления, который будет обрабатываться отдельно. Все вычисления внутри блоков будут независимы и могут выполняться параллельно. Для повышения эффективности треугольной факторизации внутренние вершины каждого из блоком можно также упорядочить с помощью метода RCM.

Более эффективными с точки зрения минимизации ширины окаймления будут следущие методы:

• Метод минимальных сепараторов, который заключается в последовательном нахождении имеющих минимальный размер разделителей (сепараторов), обеспечивающих расщепление оставшихся вершин на два независимых блока.

• Метод минимальной степени (Minimum Degree, MD). Прямое применение этого метода к матрицах большого размера затруднительно, поэтому используется приближенный метод минимальной степени (Approximate Minimum Degree, AMD).

• Метод вложенных сечений (Nested Dissection, ND). Это рекурсивный алгоритм, на каждом шаге разделяющий множество вершин на два нерависимых блока, представленые в [math]2\times2[/math] блочно-окаймленном виде.

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

Существуют и другие алгоритмы упорядочивания матриц для наиболее оптимального их распределения по процессорам. Наиболее популярными являются последовательные пакеты METIS, JOSTLE, SCOTCH, CHACO, PARTY, а также параллельные коды PARMETIS, JOSTLE, PT-SCOTCH и ZOLTAN. Многие из них являются свободно распространяемыми.

8.2 Разложение в независимых блоках

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

8.3 Разложение в сепараторах

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

8.4 Иерархические и вложенные алгоритмы

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

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

8.5 Блочный метод Якоби

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

Алтернативням является подход, когда из-за соображений увеличения ресурса параллелизма некоторые элементы отбрасываются исключительно из структурных соображений. Например, можно сначала распределить строки матрицы по процессорам, а затем перед проведением разложения просто отбросить все элементы связывающие один процессор с другим. В этом случае разложение будет проходить полностью независимо в каждом из блоков. Внутри каждого блока может проводиться любой из структурных или пороговых вариантов разложения Холецкого. Построение такого предобуславливателя называют блочным методом Якоби без перекрытия блоков (Block Jacobi, BJ). Такое предославливание является наиболее простым, применение его наиболее параллельным (полностью без обменов данными), правда сходимость может оставлять желать лучшего, почти независимо от качества разложения внутри каждого из блоков.

8.6 Адитивный метод Шварца

разложение гораздо более высокого качества по сравнению с мотодом Якоби можно получить применяя адитивный метод Шварца (Additive Schwarz, AS). Иногда этот метод называют также методом Якоби с перекрытиями. Суть его заключается в расширении структуры каждого из блоков матрицы за счет добавления нескольких слоев близлежащих строк матрицы. Треугольное разложение строится для расширенной матрицы, таким образом на каждом из процессоров происходит решение расширенной подзадачи с привлечением данных от других процессоров. После нахождения решения подзадачи на каждом из процессоров обычно происходит отбрасывание не локальных компонент решения. Такой вариант метода называют адитивный метод Шварца с ограничениями (Restricted Additive Schwarz, RAS).

Сходимость адитивного метода Шварца бывает гораздо выше чем сходимость метода Якоби, и обычно монотонно улучшается с ростом размера перекрытия. Несмотря на дополнительные вычисление и обмены, общее время решения на параллельном компьютере может быть существенно меньше.

8.7 Неполное обратное треугольное разложения

Существует и другой вариант аддитивного разложения, который кроме несколько более быстрой сходимости опирается на построение перекрытий блоков только в одну сторону ("назад", т.е. в сторону меньших номеров строк). Название этого метода блочное неполного обратного разложения Холецкого, имеющее только английскую аббревиатуру BIIC (Block Incomplete Inverse Cholesky). Позднее, вместе с рассмотрением несимметричного варианта разложения, этот метод стал называться методом неполного обратного треугольного разложения, или НОТ разложения.

Комбинация этого метода с неполным симметричным треугольным разложением второго порядка IC2 внутри каждого из блоков имеет обозначение BIIC2.

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

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

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

9.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]

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

9.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]

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

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

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

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

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

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

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

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

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