Algorithm level

Difference between revisions of "Cholesky decomposition"

From Algowiki
Jump to navigation Jump to search
[unchecked revision][checked revision]
 
(103 intermediate revisions by 5 users not shown)
Line 1: Line 1:
−
{| style="float:{{{float|right}}}; clear:right;margin-left: 1em; width:{{{width|30%}}};" cellpadding=8 cellspacing=1 border=0
+
{{algorithm
−
|-
+
| name              = Cholesky decomposition
−
|align=left width=100% style="background-color:#f3f3ff; border:1px solid"|
+
| serial_complexity = <math>O(n^3)</math>
−
'''Properties of the algorithm:'''
+
| input_data        = <math>\frac{n (n + 1)}{2}</math>
−
* Sequential complexity:<nowiki/> <math>O(n^3)</math>
+
| output_data      = <math>\frac{n (n + 1)}{2}</math>
−
* Height of the parallel form:<nowiki/> <math>O(n)</math>
+
| pf_height        = <math>O(n)</math>
−
* Width of the parallel form:<nowiki/> <math>O(n^2)</math>
+
| pf_width          = <math>O(n^2)</math>
−
* Amount of input data:<nowiki/> <math>\frac{n (n + 1)}{2}</math>
+
}}
−
* Amount of output data:<nowiki/> <math>\frac{n (n + 1)}{2}</math>
+
Main authors: [[:ru:Участник:Frolov|Alexey Frolov]].
−
|}
 
  
 
== Properties and structure of the algorithm ==
 
== Properties and structure of the algorithm ==
  
−
=== General description ===
+
=== General description of the algorithm ===
  
 +
'''The Cholesky decomposition algorithm''' was first proposed by Andre-Louis Cholesky (October 15, 1875 - August 31, 1918) at the end of the First World War shortly before he was killed in battle. He was a French military officer and mathematician. The idea of this algorithm was published in 1924 by his fellow officer<ref>Commandant Benoit, Note sur une méthode de résolution des équations normales provenant de l'application de la méthode des moindres carrés à un système d'équations linéaires en nombre inférieur à celui des inconnues (Procédé du Commandant Cholesky), Bulletin Géodésique 2 (1924), 67-77.</ref> and, later, was used by Banachiewicz in 1938<ref>Banachiewicz T. Principles d'une nouvelle technique de la méthode des moindres carrês. Bull. Intern. Acad. Polon. Sci. A., 1938,  134-135.</ref><ref>Banachiewicz T. Méthode de résoltution numérique des équations linéaires, du calcul des déterminants et des inverses et de réduction des formes quardatiques. Bull. Intern. Acad. Polon. Sci. A., 1938, 393-401.</ref>. In the Russian mathematical literature, the Cholesky decomposition is also known as the square-root method<ref>Voevodin V.V. Computational foundations of linear algebra Moscow: Nauka, 1977.</ref><ref>Voevodin V.V., Kuznetsov Yu.A. Matrices and computations, Moscow: Nauka, 1984.</ref><ref>Faddeev D.K., Faddeeva V.N. Computational methods in linear algebra. Moscow-Leningrad: Fizmatgiz, 1963.</ref> due to the square root operations used in this decomposition and not used in  Gaussian elimination.
  
−
'''The Cholesky decomposition algorithm''' was first proposed by Andre-Louis Cholesky (October 15, 1875 - August 31, 1918) at the end of the First World War shortly before he was killed in battle. He was a French military officer and mathematician. The idea of this algorithm was published in 1924 by his fellow officer and, later, was used by Banachiewicz in 1938 [7]. In the Russian mathematical literature, the Cholesky decomposition is also known as the square-root method [1-3] due to the square root operations used in this decomposition and not used in  Gaussian elimination.
+
Originally, the Cholesky decomposition was used only for dense real symmetric positive definite matrices. At present, the application of this decomposition is much wider. For example, it can also be employed for the case of Hermitian matrices. In order to increase the computing performance, its block versions are often applied.
−
 
 
−
Originally, the Cholesky decomposition was used only for dense symmetric positive definite matrices. At present, the application of this decomposition is much wider. For example, it can also be employed for the case of Hermitian matrices. In order to increase the computing performance, its block versions are often applied.
 
  
 
In the case of sparse matrices, the Cholesky decomposition is also widely used as the main stage of a direct method for solving linear systems. In order to reduce the memory requirements and the profile of the matrix, special reordering strategies are applied to minimize the number of arithmetic operations. A number of reordering strategies are used to identify the independent matrix blocks for parallel computing systems.  
 
In the case of sparse matrices, the Cholesky decomposition is also widely used as the main stage of a direct method for solving linear systems. In order to reduce the memory requirements and the profile of the matrix, special reordering strategies are applied to minimize the number of arithmetic operations. A number of reordering strategies are used to identify the independent matrix blocks for parallel computing systems.  
  
−
Various versions of the Cholesky decomposition are successfully used in iterative methods to construct preconditioners for sparse symmetric positive definite matrices. In the case of incomplete triangular decomposition, the elements of a preconditioning matrix are specified only in predetermined positions (for example, in the positions of the nonzero elements; this version is known as the IC0 decomposition). In order to construct a more accurate decomposition, a filtration of small elements is performed using a filtration threshold. The use of such a threshold allows one to obtain an accurate decomposition, but the number of nonzero elements increases. A decomposition algorithm of second-order accuracy is discussed in [11]; this algorithm retains the number of nonzero elements in the factors of the decomposition and allows one to increase the accuracy. In its parallel implementation, a special version of an additive preconditioner is applied on the basis of the second-order decomposition [12].
+
Various versions of the Cholesky decomposition are successfully used in iterative methods to construct preconditioners for sparse symmetric positive definite matrices. In the case of incomplete triangular decomposition, the elements of a preconditioning matrix are specified only in predetermined positions (for example, in the positions of the nonzero elements; this version is known as the IC0 decomposition). In order to construct a more accurate decomposition, a filtration of small elements is performed using a filtration threshold. The use of such a threshold allows one to obtain an accurate decomposition, but the number of nonzero elements increases. A decomposition algorithm of second-order accuracy is discussed in<ref>Kaporin I.E. High quality preconditioning of a general symmetric positive definite matrix based on its UTU + UTR + RTU-decomposition. Numer. Lin. Algebra Appl. (1998) Vol. 5, No. 6, 483-509.</ref>; this algorithm retains the number of nonzero elements in the factors of the decomposition and allows one to increase the accuracy. In its parallel implementation, a special version of an additive preconditioner is applied on the basis of the second-order decomposition<ref>Kaporin I.E., Kon'shin I.N. Parallel solution of symmetric positive definite systems based on decomposition into overlapping blocks. Zhurn. Vychisl. Matem. i Matem. Fiziki. (2001) Vol. 41, No. 4, 515–528.</ref>.
  
−
Here we consider the original version of the Cholesky decomposition for dense real symmetric positive definite matrices. For a number of other versions, however, the structure of this decomposition is almost the same (for example, for the complex case): the distinctions consist in the change of most real operations by the corresponding complex operations. A list of other basic versions of the Cholesky decomposition is available on the page [[The Cholesky method (symmetric triangular decomposition)]].
+
Here we consider the original version of the Cholesky decomposition for dense real symmetric positive definite matrices. For a number of other versions, however, the structure of this decomposition is almost the same (for example, for the complex case): the distinction consists in the change of the majority of real operations by the corresponding complex operations. A list of other basic versions of the Cholesky decomposition is available on the page [[Cholesky method]].
  
−
For positive definite Hermitian matrices (''symmetric matrices in the real case''), we use the decomposition <math>A = L L^*</math>, where <math>L</math> is [http://lineal.guru.ru/lineal3/texts.php?PHPSESSID=e8862lf95ipuht6hb7bd6vl1m6&&rd=3368 the lower triangular matrix], or the decomposition <math>A = U^* U</math>, where <math>U</math> is [http://lineal.guru.ru/lineal3/texts.php?PHPSESSID=e8862lf95ipuht6hb7bd6vl1m6&&rd=3368 the upper triangular matrix]. These forms of the Cholesky decomposition are equivalent in the sense of the amount of arithmetic operations and are different in the sense of data represntation. The essence of this decomposition consists in the implementation of formulas obtained uniquely for the elements of the matrix <math>L</math> form the above equality. The Cholesky decomposition is widely used due to the following features.
+
For positive definite Hermitian matrices (''symmetric matrices in the real case''), we use the decomposition <math>A = L L^*</math>, where <math>L</math> is [http://lineal.guru.ru/lineal3/texts.php?PHPSESSID=e8862lf95ipuht6hb7bd6vl1m6&&rd=3368 the lower triangular matrix], or the decomposition <math>A = U^* U</math>, where <math>U</math> is [http://lineal.guru.ru/lineal3/texts.php?PHPSESSID=e8862lf95ipuht6hb7bd6vl1m6&&rd=3368 the upper triangular matrix]. These forms of the Cholesky decomposition are equivalent in the sense of the amount of arithmetic operations and are different in the sense of data represntation. The essence of this decomposition consists in the implementation of formulas obtained uniquely for the elements of the matrix <math>L</math> from the above equality. The Cholesky decomposition is widely used due to the following features.
  
 
==== Symmetry of matrices ====
 
==== Symmetry of matrices ====
Line 34: Line 32:
  
 
The Cholesky decomposition allows one to use the so-called ''accumulation mode'' due to the fact that the significant part of computation involves ''dot product operations''. Hence, these dot products can be accumulated in double precision for additional accuracy. In this mode, the Cholesky method has the least ''[[Glossary#equivalent  perturbation | equivalent  perturbation]]''. During the process of decomposition, no growth of the matrix elements can occur, since the matrix is symmetric and positive definite.
 
The Cholesky decomposition allows one to use the so-called ''accumulation mode'' due to the fact that the significant part of computation involves ''dot product operations''. Hence, these dot products can be accumulated in double precision for additional accuracy. In this mode, the Cholesky method has the least ''[[Glossary#equivalent  perturbation | equivalent  perturbation]]''. During the process of decomposition, no growth of the matrix elements can occur, since the matrix is symmetric and positive definite.
 +
Thus, the Cholesky algorithm is unconditionally stable.
  
−
=== Mathematical description ===
+
=== Mathematical description of the algorithm ===
  
−
Исходные данные: положительно определённая симметрическая матрица <math>A</math> (элементы <math>a_{ij}</math>).
+
Input data: a symmetric positive definite matrix <math>A</math> whose elements are denoted by <math>a_{ij}</math>).
  
−
Вычисляемые данные: левая треугольная матрица <math>L</math> (элементы <math>l_{ij}</math>).
+
Output data: the lower triangular matrix <math>L</math> whose elements are denoted by <math>l_{ij}</math>).
  
−
Формулы метода:
+
The Cholesky algorithm can be represented in the form
 
:<math>
 
:<math>
 
\begin{align}
 
\begin{align}
Line 51: Line 50:
 
</math>
 
</math>
  
−
Существует также блочная версия метода, однако в данном описании разобран только точечный метод.
+
There exist block versions of this algorithm; however, here we consider only its “dot” version.
  
−
В ряде реализаций деление на диагональный элемент выполняется в два этапа: вычисление <math>\frac{1}{l_{ii}}</math> и затем умножение на него всех (видоизменённых) <math>a_{ji}</math> . Здесь мы этот вариант алгоритма не рассматриваем. Заметим только, что он имеет худшие параллельные характеристики, чем представленный.
+
In a number of implementations, the division by the diagonal element <math>l_{ii}</math> is made in the following two steps: the computation of <math>\frac{1}{l_{ii}}</math> and, then, the multiplication of the result by the modified values of <math>a_{ji}</math> . Here we do not consider this computational scheme, since this scheme has worse parallel characteristics than that given above.
  
−
=== Вычислительное ядро алгоритма ===
+
=== Computational kernel of the algorithm ===
  
−
Вычислительное ядро последовательной версии метода Холецкого можно составить из множественных (всего их <math>\frac{n (n - 1)}{2}</math>) вычислений скалярных произведений строк матрицы:
+
A computational kernel of its serial version can be composed of <math>\frac{n (n - 1)}{2}</math> dot products of the matrix rows:
  
−
:<math>\sum_{p = 1}^{i - 1} l_{ip} l_{jp}</math>
+
:<math>\sum_{p = 1}^{i - 1} l_{ip} l_{jp}.</math>
  
−
в режиме накопления или без него, в зависимости от требований задачи. В отечественных реализациях, даже в последовательных, упомянутый способ представления не используется. Дело в том, что даже в этих реализациях метода вычисление сумм типа
+
These dot products can be accumulated in double precision for additional accuracy.
 +
 
 +
In many implementations, however, the operation
  
 
:<math>a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp}</math><nowiki/>
 
:<math>a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp}</math><nowiki/>
  
−
в которых и встречаются скалярные произведения, ведутся не в порядке «вычислили скалярное произведение, а потом вычли его из элемента», а путём вычитания из элемента покомпонентных произведений, являющихся частями скалярных произведений. Поэтому следует считать вычислительным ядром метода не вычисления скалярных произведений, а вычисления выражений
+
is performed by subtracting the componentwise products as parts of dot products from <math>a_{ji}</math> instead of subtracting the entire dot product from <math>a_{ji}</math>. Hence, the following operation should be considered as a computational kernel instead of the dot product operation:
−
 
 
−
:<math>a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp}</math>
 
  
−
в ''режиме накопления'' или без него.
+
:<math>a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp}</math>.
  
−
Тем не менее, в популярных зарубежных реализациях точечного метода Холецкого, в частности, в библиотеках LINPACK и LAPACK, основанных на BLAS, используются именно вычисления скалярных произведений в виде вызова соответствующих подпрограмм BLAS (конкретно — функции SDOT). На последовательном уровне это влечёт за собой дополнительную операцию суммирования на каждый из <math>\frac{n (n + 1)}{2}</math> вычисляемый элемент матрицы <math>L</math> и некоторое замедление работы программы (о других следствиях рассказано ниже в разделе «[[#Существующие реализации алгоритма|Существующие реализации алгоритма]]»). Поэтому в данных вариантах ядром метода Холецкого будут вычисления этих скалярных произведений.
+
Here the accumulation in double precision can also be used.
  
−
=== Макроструктура алгоритма ===
+
The Cholesky algorithm is implemented in the LINPACK and LAPACK libraries on the basis of the BLAS library by using the SDOT dot product subprogram. In the serial mode, this approach involves an additional sum operation for each of the elements of the matrix <math>L</math> (the total number of these elements is equal to <math>\frac{n (n + 1)}{2}</math>) and causes a certain slowing of the algorithm performance. Other consequences are discussed in the subsection «[[#Existing implementations of the algorithm|Existing implementations of the algorithm]]»). In the BLAS-based implementations, thus, the computational kernel of the Cholesky algorithm consists of dot products.
  
−
Как записано и в [[#Вычислительное ядро алгоритма|описании ядра алгоритма]], основную часть метода составляют множественные (всего <math>\frac{n (n - 1)}{2}</math>) вычисления сумм
+
=== Macro structure of the algorithm ===
  
−
:<math>a_{ji}-\sum_{p=1}^{i-1}l_{ip} l_{jp}</math>
+
As shown above, the main part of the Cholesky algorithm consists of the following <math>\frac{n (n - 1)}{2}</math> operations (with or without accumulation):
  
−
в режиме накопления или без него.
+
:<math>a_{ji}-\sum_{p=1}^{i-1}l_{ip} l_{jp}</math>.
  
−
=== Описание схемы реализации последовательного алгоритма ===
+
=== Implementation scheme of the serial algorithm ===
  
−
Последовательность исполнения метода следующая:
+
This scheme can be represented as
  
 
1. <math>l_{11}= \sqrt{a_{11}}</math>
 
1. <math>l_{11}= \sqrt{a_{11}}</math>
  
−
2. <math>l_{j1}= \frac{a_{j1}}{l_{11}}</math> (при <math>j</math> от <math>2</math> до <math>n</math>).
+
2. <math>l_{j1}= \frac{a_{j1}}{l_{11}}</math> (for <math>j</math> from <math>2</math> to <math>n</math>).
  
−
Далее для всех <math>i</math> от <math>2</math> до <math>n</math> по нарастанию выполняются
+
For all <math>i</math> from <math>2</math> to <math>n</math>:
  
−
3. <math>l_{ii} = \sqrt{a_{ii} - \sum_{p = 1}^{i - 1} l_{ip}^2}</math> и
+
3. <math>l_{ii} = \sqrt{a_{ii} - \sum_{p = 1}^{i - 1} l_{ip}^2}</math> and
  
−
4. (кроме <math>i = n</math>): <nowiki/><math>l_{ji} = \left (a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp} \right ) / l_{ii}</math> (для всех <math>j</math> от <math>i + 1</math> до <math>n</math>).
+
4. (except for <math>i = n</math>): <nowiki/><math>l_{ji} = \left (a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp} \right ) / l_{ii}</math> (for all <math>j</math> from <math>i + 1</math> to <math>n</math>).
  
−
После этого (если <math>i < n</math>) происходит переход к шагу 3 с бо́льшим <math>i</math>.
+
If <math>i < n</math>, then <math>i = i+1</math> and go to step 3.
  
−
Особо отметим, что вычисления сумм вида <math>a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp}</math> в обеих формулах производят в режиме накопления вычитанием из <math>a_{ji}</math> произведений <math>l_{ip} l_{jp}</math> для <math>p</math> от <math>1</math> до <math>i - 1</math>, c нарастанием <math>p</math>.
+
In the above two formulas, the computation of  <math>a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp}</math> can be performed in the accumulation mode: the sum is computed in double precision and is rounded off before the subtraction.
  
−
=== Последовательная сложность алгоритма ===
+
=== Serial complexity of the algorithm ===
  
−
Для разложения матрицы порядка n методом Холецкого в последовательном (наиболее быстром) варианте требуется:
+
The following number of operations should be performed to decompose a matrix of order <math>n</math> using a serial version of the Cholesky algorithm:
 
 
−
* <math>n</math> вычислений квадратного корня,
+
* <math>n</math> square roots,
−
* <math>\frac{n(n-1)}{2}</math> делений,
+
* <math>\frac{n(n-1)}{2}</math> divisions,
−
* по <math>\frac{n^3-n}{6}</math> умножений и сложений (вычитаний) — ''основная часть алгоритма''.
+
* <math>\frac{n^3-n}{6}</math> multiplications and <math>\frac{n^3-n}{6}</math> additions (subtractions): the main amount of computational work.
−
 
−
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности (или использования функции вроде DPROD в Фортране), что ещё больше увеличивает долю умножений и сложений/вычитаний во времени, требуемом для выполнения метода Холецкого.
 
  
−
При классификации по последовательной сложности, таким образом, метод Холецкого относится к алгоритмам ''с кубической сложностью''.
+
In the accumulation mode, the multiplication and subtraction operations should be made in double precision (or by using the corresponding function, like the DPROD function in Fortran), which increases the overall computation time of the Cholesky algorithm.
  
−
=== Информационный граф ===
+
Thus, a serial version of the Cholesky algorithm is of cubic complexity.
  
−
Опишем [[глоссарий#Граф алгоритма|граф алгоритма]] как аналитически, так и в виде рисунка.
+
=== Information graph ===
  
−
Граф алгоритма состоит из трёх групп вершин, расположенных в целочисленных узлах трёх областей разной размерности.
+
The [[Glossary#algorithm graph| graph of the algorithm]]<ref>Voevodin V.V. Mathematical foundations of parallel computing. Moscow: Moscow University Press, 1991. 345 pp.</ref><ref>Voevodin V.V., Voevodin Vl.V. Parallel computing. St. Petersburg : BHV-Petersburg, St. 2002. 608 pp.</ref><ref>Frolov A.V. Design principles and description of the Sigma language. Preprint N 236. Moscow: Institute of Numerical Mathematics, 1989.</ref> consists of three groups of vertices positioned in the integer-valued nodes of three domains of different dimension.
  
−
'''Первая''' группа вершин расположена в одномерной области, соответствующая ей операция вычисляет функцию SQRT.  
+
The '''first''' group of vertices belongs to the one-dimensional domain corresponding to the SQRT operation. The only coordinate <math>i</math> of each vertex changes in the range from <math>1</math> to <math>n</math> and takes all integer values from this range.
−
Естественно введённая единственная координата каждой из вершин <math>i</math> меняется в диапазоне от <math>1</math> до <math>n</math>, принимая все целочисленные значения.
 
  
−
Аргумент этой функции
+
Argument of this operation:
 
 
−
* при <math>i = 1</math> — элемент ''входных данных'', а именно  <math>a_{11}</math>;  
+
* for <math>i = 1</math>: the argument is <math>a_{11}</math>;  
−
* при <math>i > 1</math> — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами <math>i - 1</math>, <math>i</math>, <math>i - 1</math>.
+
* for <math>i > 1</math>: the argument is the result of the operation corresponding to the vertex of the third group with coordinates <math>i - 1</math>, <math>i</math>, <math>i - 1</math>.
−
Результат срабатывания операции является ''выходным данным'' <math>l_{ii}</math>.
+
The result of the SQRT operation  is <math>l_{ii}</math>.
  
−
'''Вторая''' группа вершин расположена в одномерной области, соответствующая ей операция <math>a / b</math>.  
+
The '''second''' group of vertices belongs to the one-dimensional domain corresponding to the operation <math>a / b</math>. The coordinates of this domain are as follows:
−
Естественно введённые координаты области таковы:  
+
* the coordinate <math>i</math> changes in the range from <math>1</math> to <math>n-1</math> and takes all integer values from this range;
−
* <math>i</math> — меняется в диапазоне от <math>1</math> до <math>n-1</math>, принимая все целочисленные значения;
+
* the coordinate <math>j</math> changes in the range from <math>i+1</math> to <math>n</math> and takes all integer values from this range.
−
* <math>j</math> — меняется в диапазоне от <math>i+1</math> до <math>n</math>, принимая все целочисленные значения.
 
  
−
Аргументы операции следующие:
+
Arguments of this operation:
−
*<math>a</math>:  
+
* the numerator <math>a</math>:  
−
** при <math>i = 1</math> — элементы ''входных данных'', а именно <math>a_{j1}</math>;
+
** for <math>i = 1</math>: the element <math>a_{j1}</math>;
−
** при <math>i > 1</math> — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами <math>i - 1, j, i - 1</math>;
+
** for <math>i > 1</math>: the result of the operation corresponding to the vertex of the third group with coordinates <math>i - 1, j, i - 1</math>;
−
* <math>b</math> — результат срабатывания операции, соответствующей вершине из первой группы, с координатой <math>i</math>.
+
* the denominator <math>b</math>: the result of the operation corresponding to the vertex of the first group with coordinate <math>i</math>.
  
−
Результат срабатывания операции является ''выходным данным'' <math>l_{ji}</math>.
+
The result of this operation is <math>l_{ji}</math>.
  
−
'''Третья''' группа вершин расположена в трёхмерной области, соответствующая ей операция  <math>a - b * c</math>.  
+
The '''third''' group of vertices belongs to the three-dimensional domain corresponding to the operation <math>a - b * c</math>. The coordinates of this domain are as follows:
−
Естественно введённые координаты области таковы:  
+
* the coordinate <math>i</math> changes in the range from <math>2</math> to <math>n</math> and takes all integer values from this range;
−
* <math>i</math> — меняется в диапазоне от <math>2</math> до <math>n</math>, принимая все целочисленные значения;
+
* the coordinate <math>j</math> changes in the range from <math>i</math> to <math>n</math> and takes all integer values from this range;
−
* <math>j</math> — меняется в диапазоне от <math>i</math> до <math>n</math>, принимая все целочисленные значения;
+
* the coordinate <math>p</math> changes in the range from <math>1</math> to <math>i - 1</math> and takes all integer values from this range.
−
* <math>p</math> — меняется в диапазоне  от <math>1</math> до <math>i - 1</math>, принимая все целочисленные значения.
 
  
−
Аргументы операции следующие:
+
Arguments of this operation:
 
*<math>a</math>:
 
*<math>a</math>:
−
** при <math>p = 1</math> элемент входных данных <math>a_{ji}</math>;  
+
** for <math>p = 1</math>: the element <math>a_{ji}</math>;  
−
** при <math>p > 1</math> — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами <math>i, j, p - 1</math>;
+
** for <math>p > 1</math>: the result of the operation corresponding to the vertex of the third group with coordinates <math>i, j, p - 1</math>;
−
*<math>b</math> — результат срабатывания операции, соответствующей вершине из второй группы, с координатами <math>p, i</math>;
+
*<math>b</math>: the result of the operation corresponding to the vertex of the second group with coordinates <math>p, i</math>;
−
*<math>c</math> — результат срабатывания операции, соответствующей вершине из второй группы, с координатами <math>p, j</math>;
+
*<math>c</math>: the result of the operation corresponding to the vertex of the second group with coordinates <math>p, j</math>.
  
−
Результат срабатывания операции является ''промежуточным данным'' алгоритма.
+
The result of this operation is intermediate for the Cholesky algorithm.
  
−
Описанный граф можно посмотреть на рисунке, выполненном для случая <math>n = 4</math>. Здесь вершины первой группы обозначены жёлтым цветом и буквосочетанием sq, вершины второй — зелёным цветом и знаком деления, третьей — красным цветом и буквой f. Вершины, соответствующие операциям, производящим выходные данные алгоритма, выполнены более крупно. Дублирующие друг друга дуги даны как одна. На первом изображении показан граф алгоритма согласно классическому определению , на втором к графу алгоритма добавлены вершины , соответствующие входным данным ( обозначены синим цветом ) и выходным данным ( обозначены розовым цветом ).
+
The above graph is illustrated in Figs. 1 and 2 for <math>n = 4</math>. In these figures, the vertices of the first group are highlighted in yellow and are marked by the letters SQ; the vertices of the second group are highlighted in green and are marked by the division sign; the vertices of the third group are highlighted in red and are marked by the letter F. The vertices corresponding to the results of operations (output data) are marked by large circles. The arcs doubling one another are depicted as a single one. The representation of the graph shown in Fig.1 corresponds to the classical definition. The graph of Fig.2 contains additional vertices corresponding to input data (blue) and to output data (pink).
  
−
Можно посмотреть и [[Запись разложения Холецкого на языке Сигма|запись графа алгоритма разложения Холецкого, сделанную на языке Сигма]].
+
[[file:Cholesky full.png|thumb|center|1400px|Figure 1. The graph of the Cholesky algorithm without input and output data: SQ is the square-root operation, F is the operation a-bc, Div is division.]]
 +
[[file:Cholesky full_in_out.png|thumb|center|1400px| Figure 2. The graph of the Cholesky algorithm with input and output data: SQ is the square-root operation, F is the operation a-bc, Div is division, In and Out indicate input and output data.]]
  
−
[[file:Cholesky full.png|thumb|center|1400px|Граф алгоритма без отображения входных и выходных данных. SQ - вычисление квадратного корня, F - операция a-bc, Div - деление.]]
+
=== Parallelization resource of the algorithm ===
−
[[file:Cholesky full_in_out.png|thumb|center|1400px|Граф алгоритма с отображением входных и выходных данных. SQ - вычисление квадратного корня, F - операция a-bc, Div - деление, In - входные данные, Out - результаты.]]
 
  
−
=== Описание ресурса параллелизма алгоритма ===
+
In order to decompose a matrix of order <math>n</math>, a parallel version of the Cholesky algorithm should serially perform the following layers of operations:
 +
* <math>n</math> layers for square roots (a single square root on each layer);
 +
* <math>n - 1</math> layers for divisions (on each of the layers, the number of divisions is linear and, depending on a particular layer, varies from <math>1</math> to <math>n - 1</math>);
 +
* <math>n - 1</math> layers of multiplications and <math>n - 1</math> layers of addition/subtraction (on each of the layers, the number of these operations is quadratic and varies from <math>1</math> to <math>\frac{n^2 - n}{2}</math>).
  
−
Для разложения матрицы порядка <math>n</math> методом Холецкого в параллельном варианте требуется последовательно выполнить следующие ярусы:
+
Contrary to a serial version, in a parallel version the square-root and division operations require a significant part of overall computational time. The existence of isolated square roots on some layers of the parallel form may cause other difficulties for particular parallel computing architectures. In the case of programmable logic devices (PLD), for example, the operations of division, multiplication and addition/subtraction can be conveyorized, which is efficient in resource saving for programmable boards, whereas the isolated square roots acquire resources of such boards and force them to be out of action for a significant period of time. In the case of symmetric linear systems, the Cholesky  decomposition is preferable compared to Gaussian elimination because of the reduction in computational time by a factor of two. However, this is not true in the case of its parallel version.  
−
* <math>n</math> ярусов с вычислением квадратного корня (единичные вычисления в каждом из ярусов),
 
−
* <math>n - 1</math> ярус делений (в каждом из ярусов линейное количество делений, в зависимости от яруса — от <math>1</math> до <math>n - 1</math>),
 
−
* по <math>n - 1</math> ярусов умножений и сложений/вычитаний (в каждом из ярусов — квадратичное количество операций, от <math>1</math> до <math>\frac{n^2 - n}{2}</math>.
 
−
 
−
Таким образом, в параллельном варианте, в отличие от последовательного, вычисления квадратных корней и делений будут определять довольно значительную долю требуемого времени. При реализации на конкретных архитектурах наличие в отдельных ярусах [[глоссарий#Ярусно-параллельная форма графа алгоритма|ЯПФ]] отдельных вычислений квадратных корней может породить и другие проблемы. Например, при реализации на ПЛИСах остальные вычисления (деления и тем более умножения и сложения/вычитания) могут быть конвейеризованы, что даёт экономию и по ресурсам на программируемых платах; вычисления же квадратных корней из-за их изолированности приведут к занятию ресурсов на платах, которые будут простаивать большую часть времени. Таким образом, общая экономия в 2 раза, из-за которой метод Холецкого предпочитают в случае симметричных задач тому же методу Гаусса, в параллельном случае уже имеет место вовсе не по всем ресурсам, и главное - не по требуемому времени.
 
  
−
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности, а в параллельном варианте это означает, что практически все промежуточные вычисления для выполнения метода Холецкого в режиме накопления должны быть двойной точности. В отличие от последовательного варианта это означает увеличение требуемой памяти почти в 2 раза.
+
In addition, we should mention the fact that the accumulation mode requires multiplications and subtraction in double precision. In a parallel version, this means that almost all intermediate computations should be performed with data given in their double precision format. Contrary to a serial version, hence, this almost doubles the memory expenditure.
  
−
При классификации по высоте ЯПФ, таким образом, метод Холецкого относится к алгоритмам ''с линейной сложностью''. При классификации по ширине ЯПФ его сложность будет ''квадратичной''.
+
Thus, the Cholesky decomposition belongs to the class of algorithms of linear complexity in the sense of the height of its parallel form, whereas its complexity is quadratic in the sense of the width of its parallel form.
  
−
=== Описание входных и выходных данных ===
+
=== Input and output data of the algorithm ===
  
−
'''Входные данные''': плотная матрица <math>A</math> (элементы <math>a_{ij}</math>).
+
'''Input data''': a dense matrix <math>A</math> of order <math>n</math>; its elements are denoted by <math>a_{ij}</math>.
−
Дополнительные ограничения:
+
The additional conditions:
−
* <math>A</math> – симметрическая матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.
+
* <math>A</math> is a symmetric matrix (i.e., <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>).
−
* <math>A</math> – положительно определённая матрица, т. е. для любых ненулевых векторов <math>\vec{x}</math> выполняется <math>\vec{x}^T A \vec{x} > 0</math>.
+
* <math>A</math> is a positive definite matrix (i.e., <math>\vec{x}^T A \vec{x} > 0</math> for any nonzero vector <math>\vec{x}</math> of length <math>n</math>.
  
−
'''Объём входных данных''': <math>\frac{n (n + 1)}{2}</math> (в силу симметричности достаточно хранить только диагональ и над/поддиагональные элементы). В разных реализациях эта экономия хранения может быть выполнена разным образом. Например, в библиотеке, реализованной в НИВЦ МГУ, матрица A хранилась в одномерном массиве длины <math>\frac{n (n + 1)}{2}</math> по строкам своей нижней части.  
+
'''Amount of input data''': <math>\frac{n (n + 1)}{2}</math> ; since the matrix is symmetric, it is sufficient to store only its main diagonal and the elements above or under this diagonal. In practice, this storage saving scheme can be implemented in various ways. In the Numerical Analysis Library developed at the Moscow University Research Computing Center, for example, the lower triangle of the matrix <math>A</math> is stored row-wise in a one-dimensional array of length <math>\frac{n (n + 1)}{2}</math>.  
  
−
'''Выходные данные''': левая треугольная матрица <math>L</math> (элементы <math>l_{ij}</math>).
+
'''Output data''': the lower triangular matrix <math>L</math> whose elements are denoted by <math>l_{ij}</math>.
  
−
'''Объём выходных данных''': <math>\frac{n (n + 1)}{2}</math> (в силу треугольности достаточно хранить только ненулевые - диагональ и поддиагональные - элементы). В разных реализациях эта экономия хранения может быть выполнена разным образом. Например, в той же  библиотеке, созданной в НИВЦ МГУ, матрица <math>L</math> хранилась в одномерном массиве длины <math>\frac{n (n + 1)}{2}</math> по строкам своей нижней части.
+
'''Amount of output data''': <math>\frac{n (n + 1)}{2}</math> ; since the matrix <math>L</math> is lower triangular, it is sufficient to store only its main diagonal and the elements under this diagonal. In practice, this storage saving scheme can be implemented in various ways. In the above-mentioned library, for example, the matrix <math>L</math> is stored row-wise in a one-dimensional array of length <math>\frac{n (n + 1)}{2}</math>.
  
−
=== Свойства алгоритма ===
+
=== Properties of the algorithm ===
  
−
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является ''квадратичным'' (отношение кубической к линейной).  
+
In the case of unlimited computer resources, the ratio of the serial complexity to the parallel complexity is ''quadratic''.  
  
−
При этом вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных – всего лишь ''линейна''.
+
The computational power of the Cholesky algorithm considered as the ratio of the number of operations to the amount of input and output data is only ''linear''.
  
−
При этом алгоритм почти полностью детерминирован, это гарантируется теоремой о единственности разложения. Использование другого порядка выполнения ассоциативных операций может привести к накоплению ошибок округления, однако это влияние в используемых вариантах алгоритма не так велико, как, скажем, отказ от использования режима накопления.
+
The Cholesky is almost completely deterministic, which is ensured by the uniqueness theorem for this particular decomposition. Another order of associative operations may lead to the accumulation of round-off errors; however, the effect of this accumulation is not so large as in the case of not using the accumulation mode when computing dot products.
  
−
Дуги информационного графа, исходящие из вершин, соответствующих операциям квадратного корня и деления, образуют пучки т. н. рассылок линейной мощности (то есть степень исхода этих вершин и мощность работы с этими данными — линейная функция от порядка матрицы и координат этих вершин). При этом естественно наличие в этих пучках «длинных» дуг. Остальные дуги локальны.  
+
The information graph arcs from the vertices corresponding to the square-root and division operations can be considered as groups of data such that the function relating the multiplicity of these vertices and the number of these operations is a linear function of the matrix order and the vertex coordinates. These groups may contain «long» arcs; the remaining arcs are local.
  
−
Наиболее известной является компактная укладка графа — его проекция на треугольник матрицы, который перевычисляется укладываемыми операциями. При этом «длинные» дуги можно убрать, заменив более дальнюю пересылку комбинацией нескольких ближних (к соседям).
+
The most known is the compact packing of a graph in the form of its projection onto the matrix triangle whose elements are recomputed by the packed operations. The «long» arcs can be eliminated and replaced by a combination of short-range arcs.
  
−
[[Глоссарий#Эквивалентное возмущение|Эквивалентное возмущение]] <math>M</math> у метода Холецкого всего вдвое больше, чем возмущение <math>\delta A</math>, вносимое в матрицу при вводе чисел в компьютер:
+
The following inequality is valid for the [[Glossary#equivalent perturbation| equivalent perturbation]] <math>M</math> of the Cholesky decomposition:  
 
<math>
 
<math>
−
||M||_{E} \leq 2||\delta A||_{E}
+
||M||_{E} \leq 2||\delta A||_{E},
 
</math>
 
</math>
 +
where <math>\delta A</math> is the perturbation of the matrix <math>A</math> caused by the representation of the matrix elements in the computer memory. Such a slow growth of matrix elements during decomposition is due to the fact that the matrix is symmetric and positive definite.
  
 +
== Software implementation of the algorithm ==
  
 +
=== Implementation peculiarities of the serial algorithm ===
  
−
Это явление обусловлено положительной определённостью матрицы. Среди всех используемых разложений матриц это наименьшее из эквивалентных возмущений.
+
In its simplest version without permuting the summation, the Cholesky decomposition can be represented in Fortran as
−
 
 
−
== Программная реализация ==
 
−
 
 
−
=== Особенности реализации последовательного алгоритма ===
 
−
 
 
−
В простейшем (без перестановок суммирования) варианте метод Холецкого на Фортране можно записать так:
 
 
<source lang="fortran">
 
<source lang="fortran">
 
DO  I = 1, N
 
DO  I = 1, N
 
S = A(I,I)
 
S = A(I,I)
 
DO  IP=1, I-1
 
DO  IP=1, I-1
−
S = S - DPROD(A(I,P), A(I,IP))
+
S = S - DPROD(A(I,IP), A(I,IP))
 
END DO
 
END DO
 
A(I,I) = SQRT (S)
 
A(I,I) = SQRT (S)
Line 223: Line 212:
 
S = A(J,I)
 
S = A(J,I)
 
DO  IP=1, I-1
 
DO  IP=1, I-1
−
S = S - DPROD(A(I,P), A(J,IP))
+
S = S - DPROD(A(I,IP), A(J,IP))
 
END DO
 
END DO
 
A(J,I) = S/A(I,I)
 
A(J,I) = S/A(I,I)
Line 229: Line 218:
 
END DO
 
END DO
 
</source>
 
</source>
−
При этом для реализации режима накопления переменная <math>S</math> должна быть двойной точности.  
+
Here <math>S</math> is a double precision variable in the case of the accumulation mode.  
  
−
Для метода Холецкого существует блочная версия, которая отличается от точечной не тем, что операции над числами заменены на аналоги этих операций над блоками; её построение основано на том, что практически все циклы точечной версии имеют тип SchedDo в терминах методологии, основанной на исследовании информационного графа и, следовательно, могут быть расщеплены на составляющие. Тем не менее, обычно блочную версию метода Холецкого записывают не в виде программы с расщеплёнными и переставленными циклами, а в виде программы, подобной реализации точечного метода, в которой вместо соответствующих скалярных операций присутствуют операции над блоками.
+
A block version of the Cholesky algorithm is usually implemented in such a way that the scalar operations in its serial versions are replaced by the corresponding block-wise operations instead of using the loop unrolling and reordering techniques.
  
−
Для обеспечения локальности работы с памятью представляется более эффективной такая схема метода Холецкого (полностью эквивалентная описанной), когда исходная матрица и её разложение хранятся не в нижнем-левом, а в правом-верхнем треугольнике. Это связано с особенностью размещения массивов в Фортране и тем, что в этом случае вычисления скалярных произведений будут идти с выборкой идущих подряд элементов массива.  
+
In order to ensure the locality of memory access in the Cholesky algorithm, in its Fortran implementation the original matrix and its decomposition are stored in the upper triangle instead of the lower triangle. The efficiency of such a version can be explained by the fact that Fortran stores matrices by columns and, hence, the computer programs in which the inner loops go up or down a column generate serial access to memory, contrary to the non-serial access when the inner loop goes across a row. This column orientation provides a significant improvement on computers with paging and cache memory.
  
−
Есть и другой вариант точечной схемы: использовать вычисляемые элементы матрицы <math>L</math> в качестве аргументов непосредственно «сразу после» их вычисления. Такая программа будет выглядеть так:
+
There exists the following dot version of the Cholesky decomposition: the elements of the matrix <math>L</math> are used as arguments in subsequent operations as soon as they are computed. This version can be illustrated as follows:
 
<source lang="fortran">
 
<source lang="fortran">
 
DO  I = 1, N
 
DO  I = 1, N
Line 249: Line 238:
 
END DO
 
END DO
 
</source>
 
</source>
−
Как видно, в этом варианте для реализации режима накопления одинарной точности мы должны будем объявить двойную точность для массива, хранящего исходные данные и результат. Подчеркнём, что [[глоссарий#Граф алгоритма|граф алгоритма]] обеих схем - один и тот же (из п.1.7), если не считать изменением замену умножения на функцию DPROD!
+
As can be seen from the above program fragment, the array to store the original matrix and the output data should be declared as double precision for the accumulation mode. Note that the [[Glossary#graph of the algorithm| graph of the algorithm]] for this fragment and for the previous one is almost the same (the only distinction is that the DPROD function is used instead of multiplications).
−
 
 
−
=== Описание локальности данных и вычислений ===
 
−
 
 
−
==== Описание локальности алгоритма ====
 
−
 
 
−
==== Описание локальности реализации алгоритма ====
 
−
 
 
−
===== Описание структуры обращений в память и качественная оценка локальности =====
 
−
 
 
−
[[file:Cholesky_locality1.jpg|thumb|center|700px|Рисунок 1. Реализация метода Холецкого. Общий профиль обращений в память]]
 
−
 
 
−
На рисунке 1 представлен профиль обращений в память для реализации метода Холецкого. В программе задействован только 1 массив, поэтому в данном случае обращения в профиле происходят только к элементам этого массива. Программа состоит из одного основного этапа, который, в свою очередь, состоит из последовательности подобных итераций. Пример одной итерации выделен зеленым цветом.
 
−
 
 
−
Видно, что на каждой <math>i</math>-й итерации используются все адреса, кроме первых <math>k_i</math>, при этом с ростом <math>i</math> увеличивается значение <math>k_i</math>. Также можно заметить, что число обращений в память на каждой итерации растет примерно до середины работы программы, после чего уменьшается вплоть до завершения работы. Это позволяет говорить о том, что данные в программе используются неравномерно, при этом многие итерации, особенно в начале выполнения программы, задействуют большой объем данных, что приводит к ухудшению локальности.
 
−
 
 
−
Однако в данном случае основным фактором, влияющим на локальность работы с памятью, является строение итерации. Рассмотрим фрагмент профиля, соответствующий нескольким первым итерациям.
 
−
 
 
−
[[file:Cholesky_locality2.jpg|thumb|center|700px|Рисунок 2. Реализация метода Холецкого. Фрагмент профиля (несколько первых итераций)]]
 
−
 
 
−
Исходя из рисунка 2 видно, что каждая итерация состоит из двух различных фрагментов. Фрагмент 1 – последовательный перебор (с некоторым шагом) всех адресов, начиная с некоторого начального. При этом к каждому элементу происходит мало обращений. Такой фрагмент обладает достаточно неплохой пространственной локальностью, так как шаг по памяти между соседними обращениями невелик, но плохой временно́й локальностью, поскольку данные редко используются повторно.
 
−
 
 
−
Фрагмент 2 устроен гораздо лучше с точки зрения локальности. В рамках этого фрагмента выполняется большое число обращений подряд к одним и тем же данным, что обеспечивает гораздо более высокую степень как пространственной, так и временно́й локальности по сравнению с фрагментом 1.
 
−
 
 
−
После рассмотрения фрагмента профиля на рис. 2 можно оценить общую локальность двух фрагментов на каждой итерации. Однако стоит рассмотреть более подробно, как устроен каждый из фрагментов.
 
−
 
 
−
[[file:Cholesky_locality3.jpg|thumb|center|700px|Рисунок 3. Реализация метода Холецкого. Фрагмент профиля (часть одной итерации)]]
 
−
 
 
−
Рис. 3, на котором представлена часть одной итерации общего профиля,  позволяет отметить достаточно интересный факт: строение каждого из фрагментов на самом деле заметно сложнее, чем это выглядит на рис. 2. В частности, каждый шаг фрагмента 1 состоит из нескольких обращений к соседним элементам, причем выполняется не последовательный перебор. Также можно увидеть, что фрагмент 2 на самом деле в свою очередь состоит из повторяющихся итераций, при этом видно, что каждый шаг фрагмента 1 соответствует одной итерации фрагмента 2 (выделено зеленым на рис. 3). Это лишний раз говорит о том, что для точного понимания локальной структуры профиля необходимо его рассмотреть на уровне отдельных обращений.
 
−
 
 
−
Стоит отметить, что выводы на основе рис. 3 просто дополняют общее представлении о строении профиля обращений; сделанные на основе рис. 2 выводы относительно общей локальности двух фрагментов остаются верны.
 
−
 
 
−
===== Количественная оценка локальности =====
 
  
−
Основной фрагмент реализации, на основе которого были получены количественные оценки, приведен [http://git.algowiki-project.org/Voevodin/locality/blob/master/benchmarks/holecky/holecky.h здесь] (функция Kernel). Условия запуска описаны [http://git.algowiki-project.org/Voevodin/locality/blob/master/README.md здесь].
+
=== Possible methods and considerations for parallel implementation of the algorithm  ===
  
−
Первая оценка выполняется на основе характеристики daps, которая оценивает число выполненных обращений (чтений и записей) в память в секунду. Данная характеристика является аналогом оценки flops применительно к работе с памятью и является в большей степени оценкой производительности взаимодействия с памятью, чем оценкой локальности. Однако она служит хорошим источником информации, в том числе для сравнения с результатами по следующей характеристике cvg.
+
An analysis of the algorithm’s structure allows one to conclude that a large variety of its parallel implementations can be proposed. In the second version (see «[[#Implementation peculiarities of the serial algorithm|Implementation peculiarities of the serial algorithm]]»), for example, all the inner loops are parallel, whereas, in the first version, only the inner loop over <math>J</math> is parallel. Nevertheless, a simple parallelization technique causes a large number of data transfer between the processors at each step of the outer loop; this number is almost comparable with the number of arithmetic operations. Hence, it is reasonable to partition the computations into blocks with the corresponding partitioning of the data arrays before the allocation of operations and data between the processors of the computing system in use.
  
−
[[file:Cholesky_locality4.jpg|thumb|center|700px|Рисунок 4. Сравнение значений оценки daps]]
+
A good-enough decision depends on the features of a particular computer system. If the nodes of a multiprocessor computer are equipped with conveyors, it is reasonable to compute several dot products at once in parallel. Such a possibility exists in the case of programmable logic devices; in this case, however, the arithmetic speed is limited by a slow serial square-root operation. In principle, it is possible to apply the so-called «skew» parallelism; however, this approach is not used in practice because of complexity in the control structure of the algorithm’s implementation.
  
−
На рисунке 4 приведены значения daps для реализаций распространенных алгоритмов, отсортированные по возрастанию (чем больше daps, тем в общем случае выше производительность). Можно увидеть, что реализация метода Холецкого характеризуется достаточно высокой скоростью взаимодействия с памятью, однако ниже, чем, например, у теста Линпак или реализации метода Якоби.
+
The dot version of the Cholesky algorithm is implemented in the basic libraries produced in Russia and in the libraries produced in western countries (for example, in LINPACK, LAPACK, ScaLAPACK, and others).
  
−
Вторая характеристика – cvg – предназначена для получения более машинно-независимой оценки локальности. Она определяет, насколько часто в программе необходимо подтягивать данные в кэш-память. Соответственно, чем меньше значение cvg, тем реже это нужно делать, тем лучше локальность.  
+
In the Russian libraries, as a rule, the accumulation mode is implemented to reduce the effect of round-off errors. In order to save computer memory, a packed representation of the matrices <math>A</math> and <math>L</math> is also used in the form of one-dimensional arrays.
  
−
[[file:Cholesky_locality5.jpg|thumb|center|700px|Рисунок 5. Сравнение значений оценки cvg]]
+
In the modern western libraries, the dot Cholesky implementations are based on its LINPACK implementation utilizing the BLAS library. The dot product is computed in BLAS with no accumulation mode, but this mode can easily be implemented by minor modifications of the SDOT subprogram included in the BLAS library.
  
−
На рисунке 5 приведены значения cvg для того же набора реализаций, отсортированные по убыванию (чем меньше cvg, тем в общем случае выше локальность). Можно увидеть, что, согласно данной оценке, реализация метода Холецкого оказалась ниже в списке по сравнению с оценкой daps.
+
It is interesting to note that the original Cholesky decomposition is available in the largest western libraries, whereas the faster <math>LU</math> decomposition algorithm without square-root operations is used only in special cases (for example, for tridiagonal matrices) when the number of diagonal elements is comparable with the number of off-diagonal elements.
  
−
<!--
+
=== Run results ===
−
===== Анализ на основе теста Apex-Map =====
 
−
--->
 
  
−
=== Возможные способы и особенности реализации параллельного алгоритма ===
+
=== Conclusions for different classes of computer architecture ===
  
−
Как нетрудно видеть по структуре графа алгоритма, вариантов распараллеливания алгоритма довольно много. Например, во втором варианте (см. раздел «[[#Особенности реализации последовательного алгоритма|Особенности реализации последовательного алгоритма]]») все внутренние циклы параллельны, в первом — параллелен цикл по <math>J</math>. Тем не менее, простое распараллеливание таким способом «в лоб» вызовет такое количество пересылок между процессорами с каждым шагом по внешнему циклу, которое почти сопоставимо с количеством арифметических операций. Поэтому перед размещением операций и данных между процессорами вычислительной системы предпочтительно разбиение всего пространства вычислений на блоки, с сопутствующим разбиением обрабатываемого массива.
+
Analyzing the performance of the ScaLAPACK library on supercomputers, we can conclude that, for large values of the matrix order <math>n</math>, the data exchanges reduce the execution efficiency, but to a smaller extent than the nonoptimality in the organization of computations within a single node. At the first stages, hence, it is necessary to optimize not a block algorithm but the subroutines used on individual processors, such as the dot Cholesky decomposition, matrix multiplications, etc. A number of possible directions of such an optimization are discussed [[#Existing implementations of the Cholesky algorithm|below]].
  
−
Многое зависит от конкретного типа вычислительной системы. Присутствие конвейеров на узлах многопроцессорной системы делает рентабельным параллельное вычисление нескольких скалярных произведений сразу. Подобная возможность есть и на программировании ПЛИСов, но там быстродействие будет ограничено медленным последовательным выполнением операции извлечения квадратного корня.  
+
Generally speaking, the efficiency of the Cholesky algorithm cannot be high for parallel computer architectures. This fact can explained by the following property of its information structure: the square-root operations are the bottleneck of the Cholesky algorithm in comparison with the division operations or with the <math>a - bc</math> operations, since these operations can easily be conveyorized or distributed among the nodes. In order to enhance the performance efficiency on parallel computers, hence, it is reasonable to use not the original Cholesky algorithm but its well-known modification without square-root operations; see [[A compact scheme of Gaussian decomposition (a version for symmetric matrices)|the matrix decomposition in the form <math>L D L^T</math>]]<ref>Krishnamoorthy A., Menon D. Matrix Inversion Using Cholesky Decomposition. 2013. eprint arXiv:1111.4144</ref>.
−
 
−
В принципе, возможно и использование т. н. «скошенного» параллелизма. Однако его на практике никто не использует, из-за усложнения управляющей структуры программы.
 
  
−
=== Масштабируемость алгоритма и его реализации ===
+
== References ==
−
 
 
−
==== Описание масштабируемости алгоритма ====
 
−
 
 
−
==== Описание масштабируемости реализации алгоритма ====
 
−
Проведём исследование масштабируемости параллельной реализации разложения Холецкого согласно [[Scalability methodology|методике]]. Исследование проводилось на суперкомпьютере "Ломоносов" [http://parallel.ru/cluster Суперкомпьютерного комплекса Московского университета].
 
−
 
 
−
Набор и границы значений изменяемых [[Глоссарий#Параметры запуска|параметров запуска]] реализации алгоритма:
 
−
 
 
−
* число процессоров [4 : 256] с шагом 4;
 
−
* размер матрицы [1024 : 5120].
 
−
 
 
−
В результате проведённых экспериментов был получен следующий диапазон [[Глоссарий#Эффективность реализации|эффективности реализации]] алгоритма:
 
−
 
 
−
* минимальная эффективность реализации 0,11%;
 
−
* максимальная эффективность реализации 2,65%.
 
−
 
 
−
На следующих рисунках приведены графики [[Глоссарий#Производительность|производительности]] и эффективности выбранной реализации разложения Холецкого в зависимости от изменяемых параметров запуска.
 
−
 
 
−
[[file:Масштабируемость Параллельной реализации метода Холецкого Производительность.png|thumb|center|700px|Рисунок 1. Параллельная реализация метода Холецкого. Изменение производительности в зависимости от числа процессоров и размера матрицы.]]
 
−
[[file:Холецкий масштабируемость эффективность.png|thumb|center|700px|Рисунок 2. Параллельная реализация метода Холецкого. Изменение производительности в зависимости от числа процессоров и размера матрицы.]]
 
−
 
 
−
Построим оценки масштабируемости выбранной реализации разложения Холецкого:
 
−
* По числу процессов: -0,000593. При увеличении числа процессов эффективность на рассмотренной области изменений параметров запуска уменьшается, однако в целом уменьшение не очень быстрое. Малая интенсивность изменения объясняется крайне низкой общей эффективностью работы приложения с максимумом в 2,65%, и значение эффективности на рассмотренной области значений быстро доходит до десятых долей процента. Это свидетельствует о том, что на большей части области значений нет интенсивного снижения эффективности. Это объясняется также тем, что с ростом [[Глоссарий#Вычислительная сложность|вычислительной сложности]] падение эффективности становится не таким быстрым. Уменьшение эффективности на рассмотренной области работы параллельной программы объясняется быстрым ростом накладных расходов на организацию параллельного выполнения. С ростом вычислительной сложности задачи эффективность снижается так же быстро, но при больших значениях числа процессов. Это подтверждает предположение о том, что накладные расходы начинают сильно превалировать над вычислениями.
 
−
* По размеру задачи: 0,06017. При увеличении размера задачи эффективность возрастает. Эффективность возрастает тем быстрее, чем большее число процессов используется для выполнения. Это подтверждает предположение о том, что размер задачи сильно влияет на эффективность выполнения приложения. Оценка показывает, что с ростом размера задачи эффективность на рассмотренной области значений параметров запуска сильно увеличивается. Также, учитывая разницу максимальной и минимальной эффективности в 2,5%, можно сделать вывод, что рост эффективности при увеличении размера задачи наблюдается на большей части рассмотренной области значений.
 
−
* По двум направлениям: 0,000403. При рассмотрении увеличения как вычислительной сложности, так и числа процессов на всей рассмотренной области значений эффективность увеличивается, однако скорость увеличения эффективности небольшая. В совокупности с тем фактом, что разница между максимальной и минимальной эффективностью на рассмотренной области значений параметров небольшая, эффективность с увеличением масштабов возрастает, но медленно и с небольшими перепадами.
 
−
 
 
−
 
 
−
[http://git.algowiki-project.org/Teplov/Scalability/tree/master/cholesky-decomposition-master Исследованная параллельная реализация на языке C]
 
−
 
 
−
=== Динамические характеристики и эффективность реализации алгоритма ===
 
−
 
 
−
Для существующих параллельных реализаций характерно отнесение всего ресурса параллелизма на блочный уровень. Относительно низкая эффективность работы связана с проблемами внутри одного узла, следующим фактором является неоптимальное соотношение между «арифметикой» и обменами. Видно, что при некотором (довольно большом) оптимальном размере блока обмены влияют не так уж сильно.
 
−
 
 
−
Для проведения экспериментов использовалась реализация разложения Холецкого, представленная в пакете SCALAPACK библиотеки Intel MKL (метод pdpotrf).  Все результаты получены на суперкомпьютере «Ломоносов». Использовались процессоры Intel Xeon X5570 с пиковой производительностью в 94 Гфлопс, а также компилятор Intel с опцией –O2.
 
−
 
 
−
[[file:Cholesky_efficiency1.jpg|thumb|center|700px]]
 
−
 
 
−
На рисунке показана эффективность реализации разложения Холецкого (случай использования нижних треугольников матриц) для разного числа процессов и разной размерности матрицы (10000, 20000 и 50000 элементов). Видно, что общая эффективность достаточно невысока – даже при малом числе процессов менее 10 %. Однако при всех размерностях матрицы эффективность снижается очень медленно – в самом худшем случае, при N = 10000, при увеличении числа процессов с 16 до 900 (в 56 раз) эффективность падает с 7 % до 0,8 % (всего в 9 раз). При N = 50000 эффективность уменьшается еще медленнее.
 
−
 
 
−
Также стоит отметить небольшое суперлинейное ускорение, полученное на 4-х процессах для N = 10000 и N = 20000 (для N = 50000 провести эксперименты на таком малом числе процессов не удалось). Помимо более эффективной работы с кэш-памятью, оно, видимо, обусловлено тем, что эти 4 процесса помещаются на ядрах одного процессора, что позволяет использовать общую память, и, соответственно, не приводит к появлению накладных расходов, связанных с пересылкой данных.
 
−
 
 
−
=== Выводы для классов архитектур ===
 
−
 
 
−
Как видно по показателям SCALAPACK на суперкомпьютерах, обмены при большом n хоть и уменьшают эффективность расчётов, но слабее, чем неоптимальность организации расчётов внутри одного узла. Поэтому, видимо, следует сначала оптимизировать не блочный алгоритм, а подпрограммы, используемые на отдельных процессорах: точечный метод Холецкого, перемножения матриц и др. подпрограммы. [[#Существующие реализации алгоритма|Ниже]] содержится информация о возможном направлении такой оптимизации.
 
−
 
 
−
В отношении же архитектуры типа ПЛИС вполне показателен тот момент, что разработчики — наполнители библиотек для ПЛИСов пока что не докладывают об успешных и эффективных реализациях точечного метода Холецкого на ПЛИСах. Это связано со следующим свойством информационной структуры алгоритма: если операции деления или вычисления выражений <math>a - bc</math> являются не только массовыми, но и параллельными, и потому их вычисления сравнительно легко выстраивать в конвейеры, то операции извлечения квадратных корней являются узким местом алгоритма - отведённое на эту операцию оборудование неизбежно будет простаивать большую часть времени. Поэтому для эффективной реализации на ПЛИСах алгоритмов, столь же хороших по вычислительным характеристикам, как и метод квадратного корня, следует использовать не метод Холецкого, а его давно известную модификацию без извлечения квадратных корней — [[Компактная схема разложения Гаусса (вариант для симметричных матриц)|разложение матрицы в произведение <math>L D L^T</math>]].
 
−
 
 
−
=== Существующие реализации алгоритма ===
 
−
 
 
−
Точечный метод Холецкого реализован как в основных библиотеках программ отечественных организаций, так и в западных пакетах LINPACK, LAPACK, SCALAPACK и др.
 
−
 
−
При этом в отечественных реализациях, как правило, выполнены стандартные требования к методу с точки зрения ошибок округления, то есть, реализован режим накопления, и обычно нет лишних операций. Правда, анахронизмом в наше время выглядит то, что ряд старых отечественных реализаций использует для экономии памяти упаковку матриц <math>A</math> и <math>L</math> в одномерный массив. При реальных вычислениях на современных вычислительных системах данная упаковка только создаёт дополнительные накладные расходы. Однако отечественные реализации, не использующие такую упаковку, вполне отвечают требованиям современности в отношении вычислительной точности метода. 
 
−
 
−
Реализация точечного метода Холецкого в современных западных пакетах страдает другими недостатками, вытекающими из того, что все эти реализации, по сути, происходят из одной и той же реализации метода в LINPACK, а та использует пакет BLAS. Основным их недостатком является даже не наличие лишних операций, о котором уже написано в разделе «[[#Вычислительное ядро алгоритма|Вычислительное ядро алгоритма]]», ибо их количество всё же невелико по сравнению с главной частью, а то, что в BLAS скалярное произведение реализовано без режима накопления. Это перечёркивает имеющиеся для метода Холецкого хорошие оценки влияния ошибок округления, поскольку они выведены как раз в предположении использования режима накопления при вычислении скалярных произведений. Поэтому тот, кто использует реализации метода Холецкого из LINPACK, LAPACK, SCALAPACK и т. п. пакетов, серьёзно рискует не получить требуемую вычислительную точность, либо ему придётся для получения хороших оценок ''одинарной точности'' использовать подпрограммы ''двойной''.
 
  
−
В крупнейших библиотеках алгоритмов до сих пор предлагается именно разложение Холецкого, а более быстрый алгоритм LU-разложения без извлечения квадратных корней используется только в особых случаях (например, для трёхдиагональных матриц), в которых количество диагональных элементов уже сравнимо с количеством внедиагональных.
+
<references />
  
−
== Литература ==
+
[[Ru:Разложение Холецкого (метод квадратного корня)]]
  
−
# Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977.
+
[[Category:Finished articles]]
−
# Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. М.: Наука, 1984.
 
−
# Фаддеев Д.К., Фаддева В.Н. Вычислительные основы линейной алгебры. М.-Л.: Физматгиз, 1963.
 
−
# Воеводин В.В.  Математические основы параллельных вычислений// М.: Изд. Моск. ун-та, 1991. 345 с.
 
−
# Фролов А.В.. Принципы построения и описание языка Сигма. Препринт ОВМ АН N 236. М.: ОВМ АН СССР, 1989.
 
−
# Commandant Benoit, Note sur une méthode de résolution des équations normales provenant de l'application de la méthode des moindres carrés à un système d'équations linéaires en nombre inférieur à celui des inconnues (Procédé du Commandant Cholesky), Bulletin Géodésique 2 (1924), 67-77.
 
−
# Banachiewicz T. Principles d'une nouvelle technique de la méthode des moindres carrês. Bull. Intern. Acad. Polon. Sci. A., 1938,  134-135.
 
−
# Banachiewicz T. Méthode de résoltution numérique des équations linéaires, du calcul des déterminants et des inverses et de réduction des formes quardatiques. Bull. Intern. Acad. Polon. Sci. A., 1938, 393-401.
 
−
# Krishnamoorthy A., Menon D. Matrix Inversion Using Cholesky Decomposition. 2013. eprint arXiv:1111.4144
 
−
# Открытая энциклопедия свойств алгоритмов. URL: http://algowiki-project.org
 
−
# Kaporin I.E. High quality preconditioning of a general symmetric positive definite matrix based on its UTU + UTR + RTU-decomposition. Numer. Lin. Algebra Appl. (1998) Vol. 5, No. 6, 483-509.
 
−
# Капорин И.Е., Коньшин И.Н. Параллельное решение симметричных положительно-определенных систем на основе перекрывающегося разбиения на блоки. Ж. вычисл. матем. и матем. физ.,  2001, Т, 41, N. 4, C. 515–528.
 
−
# Воеводин Вл., Жуматий С., Соболев С., Антонов А., Брызгалов П., Никитенко Д., Стефанов К., Воеводин Вад. Практика суперкомпьютера «Ломоносов» // Открытые системы, 2012, N 7, С. 36-39.
 
−
# Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. – СПб.: БХВ - Петербург, 2002. – 608 с.
 
−
# Воеводин Вад. В. Визуализация и анализ профиля обращений в память // Вестник Южно-Уральского государственного университета. Серия Математическое моделирование и про-граммирование. — 2011. — Т. 17, № 234. — С. 76–84.
 
−
# Воеводин Вл. В., Воеводин Вад. В. Спасительная локальность суперкомпьютеров // Откры-тые системы. — 2013. — № 9. — С. 12–15.
 

Latest revision as of 18:37, 2 July 2023


Cholesky decomposition
Sequential algorithm
Serial complexity [math]O(n^3)[/math]
Input data [math]\frac{n (n + 1)}{2}[/math]
Output data [math]\frac{n (n + 1)}{2}[/math]
Parallel algorithm
Parallel form height [math]O(n)[/math]
Parallel form width [math]O(n^2)[/math]

Main authors: Alexey Frolov.

1 Properties and structure of the algorithm

1.1 General description of the algorithm

The Cholesky decomposition algorithm was first proposed by Andre-Louis Cholesky (October 15, 1875 - August 31, 1918) at the end of the First World War shortly before he was killed in battle. He was a French military officer and mathematician. The idea of this algorithm was published in 1924 by his fellow officer[1] and, later, was used by Banachiewicz in 1938[2][3]. In the Russian mathematical literature, the Cholesky decomposition is also known as the square-root method[4][5][6] due to the square root operations used in this decomposition and not used in Gaussian elimination.

Originally, the Cholesky decomposition was used only for dense real symmetric positive definite matrices. At present, the application of this decomposition is much wider. For example, it can also be employed for the case of Hermitian matrices. In order to increase the computing performance, its block versions are often applied.

In the case of sparse matrices, the Cholesky decomposition is also widely used as the main stage of a direct method for solving linear systems. In order to reduce the memory requirements and the profile of the matrix, special reordering strategies are applied to minimize the number of arithmetic operations. A number of reordering strategies are used to identify the independent matrix blocks for parallel computing systems.

Various versions of the Cholesky decomposition are successfully used in iterative methods to construct preconditioners for sparse symmetric positive definite matrices. In the case of incomplete triangular decomposition, the elements of a preconditioning matrix are specified only in predetermined positions (for example, in the positions of the nonzero elements; this version is known as the IC0 decomposition). In order to construct a more accurate decomposition, a filtration of small elements is performed using a filtration threshold. The use of such a threshold allows one to obtain an accurate decomposition, but the number of nonzero elements increases. A decomposition algorithm of second-order accuracy is discussed in[7]; this algorithm retains the number of nonzero elements in the factors of the decomposition and allows one to increase the accuracy. In its parallel implementation, a special version of an additive preconditioner is applied on the basis of the second-order decomposition[8].

Here we consider the original version of the Cholesky decomposition for dense real symmetric positive definite matrices. For a number of other versions, however, the structure of this decomposition is almost the same (for example, for the complex case): the distinction consists in the change of the majority of real operations by the corresponding complex operations. A list of other basic versions of the Cholesky decomposition is available on the page Cholesky method.

For positive definite Hermitian matrices (symmetric matrices in the real case), we use the decomposition [math]A = L L^*[/math], where [math]L[/math] is the lower triangular matrix, or the decomposition [math]A = U^* U[/math], where [math]U[/math] is the upper triangular matrix. These forms of the Cholesky decomposition are equivalent in the sense of the amount of arithmetic operations and are different in the sense of data represntation. The essence of this decomposition consists in the implementation of formulas obtained uniquely for the elements of the matrix [math]L[/math] from the above equality. The Cholesky decomposition is widely used due to the following features.

1.1.1 Symmetry of matrices

The symmetry of a matrix allows one to store in computer memory slightly more than half the number of its elements and to reduce the number of operations by a factor of two compared to Gaussian elimination. Note that the LU-decomposition does not require the square-root operations when using the property of symmetry and, hence, is somewhat faster than the Cholesky decomposition, but requires to store the entire matrix.

1.1.2 Accumulation mode

The Cholesky decomposition allows one to use the so-called accumulation mode due to the fact that the significant part of computation involves dot product operations. Hence, these dot products can be accumulated in double precision for additional accuracy. In this mode, the Cholesky method has the least equivalent perturbation. During the process of decomposition, no growth of the matrix elements can occur, since the matrix is symmetric and positive definite. Thus, the Cholesky algorithm is unconditionally stable.

1.2 Mathematical description of the algorithm

Input data: a symmetric positive definite matrix [math]A[/math] whose elements are denoted by [math]a_{ij}[/math]).

Output data: the lower triangular matrix [math]L[/math] whose elements are denoted by [math]l_{ij}[/math]).

The Cholesky algorithm can be represented in the form

[math] \begin{align} l_{11} & = \sqrt{a_{11}}, \\ l_{j1} & = \frac{a_{j1}}{l_{11}}, \quad j \in [2, n], \\ l_{ii} & = \sqrt{a_{ii} - \sum_{p = 1}^{i - 1} l_{ip}^2}, \quad i \in [2, n], \\ l_{ji} & = \left (a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp} \right ) / l_{ii}, \quad i \in [2, n - 1], j \in [i + 1, n]. \end{align} [/math]

There exist block versions of this algorithm; however, here we consider only its “dot” version.

In a number of implementations, the division by the diagonal element [math]l_{ii}[/math] is made in the following two steps: the computation of [math]\frac{1}{l_{ii}}[/math] and, then, the multiplication of the result by the modified values of [math]a_{ji}[/math] . Here we do not consider this computational scheme, since this scheme has worse parallel characteristics than that given above.

1.3 Computational kernel of the algorithm

A computational kernel of its serial version can be composed of [math]\frac{n (n - 1)}{2}[/math] dot products of the matrix rows:

[math]\sum_{p = 1}^{i - 1} l_{ip} l_{jp}.[/math]

These dot products can be accumulated in double precision for additional accuracy.

In many implementations, however, the operation

[math]a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp}[/math]

is performed by subtracting the componentwise products as parts of dot products from [math]a_{ji}[/math] instead of subtracting the entire dot product from [math]a_{ji}[/math]. Hence, the following operation should be considered as a computational kernel instead of the dot product operation:

[math]a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp}[/math].

Here the accumulation in double precision can also be used.

The Cholesky algorithm is implemented in the LINPACK and LAPACK libraries on the basis of the BLAS library by using the SDOT dot product subprogram. In the serial mode, this approach involves an additional sum operation for each of the elements of the matrix [math]L[/math] (the total number of these elements is equal to [math]\frac{n (n + 1)}{2}[/math]) and causes a certain slowing of the algorithm performance. Other consequences are discussed in the subsection «Existing implementations of the algorithm»). In the BLAS-based implementations, thus, the computational kernel of the Cholesky algorithm consists of dot products.

1.4 Macro structure of the algorithm

As shown above, the main part of the Cholesky algorithm consists of the following [math]\frac{n (n - 1)}{2}[/math] operations (with or without accumulation):

[math]a_{ji}-\sum_{p=1}^{i-1}l_{ip} l_{jp}[/math].

1.5 Implementation scheme of the serial algorithm

This scheme can be represented as

1. [math]l_{11}= \sqrt{a_{11}}[/math]

2. [math]l_{j1}= \frac{a_{j1}}{l_{11}}[/math] (for [math]j[/math] from [math]2[/math] to [math]n[/math]).

For all [math]i[/math] from [math]2[/math] to [math]n[/math]:

3. [math]l_{ii} = \sqrt{a_{ii} - \sum_{p = 1}^{i - 1} l_{ip}^2}[/math] and

4. (except for [math]i = n[/math]): [math]l_{ji} = \left (a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp} \right ) / l_{ii}[/math] (for all [math]j[/math] from [math]i + 1[/math] to [math]n[/math]).

If [math]i \lt n[/math], then [math]i = i+1[/math] and go to step 3.

In the above two formulas, the computation of [math]a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp}[/math] can be performed in the accumulation mode: the sum is computed in double precision and is rounded off before the subtraction.

1.6 Serial complexity of the algorithm

The following number of operations should be performed to decompose a matrix of order [math]n[/math] using a serial version of the Cholesky algorithm:

  • [math]n[/math] square roots,
  • [math]\frac{n(n-1)}{2}[/math] divisions,
  • [math]\frac{n^3-n}{6}[/math] multiplications and [math]\frac{n^3-n}{6}[/math] additions (subtractions): the main amount of computational work.

In the accumulation mode, the multiplication and subtraction operations should be made in double precision (or by using the corresponding function, like the DPROD function in Fortran), which increases the overall computation time of the Cholesky algorithm.

Thus, a serial version of the Cholesky algorithm is of cubic complexity.

1.7 Information graph

The graph of the algorithm[9][10][11] consists of three groups of vertices positioned in the integer-valued nodes of three domains of different dimension.

The first group of vertices belongs to the one-dimensional domain corresponding to the SQRT operation. The only coordinate [math]i[/math] of each vertex changes in the range from [math]1[/math] to [math]n[/math] and takes all integer values from this range.

Argument of this operation:

  • for [math]i = 1[/math]: the argument is [math]a_{11}[/math];
  • for [math]i \gt 1[/math]: the argument is the result of the operation corresponding to the vertex of the third group with coordinates [math]i - 1[/math], [math]i[/math], [math]i - 1[/math].

The result of the SQRT operation is [math]l_{ii}[/math].

The second group of vertices belongs to the one-dimensional domain corresponding to the operation [math]a / b[/math]. The coordinates of this domain are as follows:

  • the coordinate [math]i[/math] changes in the range from [math]1[/math] to [math]n-1[/math] and takes all integer values from this range;
  • the coordinate [math]j[/math] changes in the range from [math]i+1[/math] to [math]n[/math] and takes all integer values from this range.

Arguments of this operation:

  • the numerator [math]a[/math]:
    • for [math]i = 1[/math]: the element [math]a_{j1}[/math];
    • for [math]i \gt 1[/math]: the result of the operation corresponding to the vertex of the third group with coordinates [math]i - 1, j, i - 1[/math];
  • the denominator [math]b[/math]: the result of the operation corresponding to the vertex of the first group with coordinate [math]i[/math].

The result of this operation is [math]l_{ji}[/math].

The third group of vertices belongs to the three-dimensional domain corresponding to the operation [math]a - b * c[/math]. The coordinates of this domain are as follows:

  • the coordinate [math]i[/math] changes in the range from [math]2[/math] to [math]n[/math] and takes all integer values from this range;
  • the coordinate [math]j[/math] changes in the range from [math]i[/math] to [math]n[/math] and takes all integer values from this range;
  • the coordinate [math]p[/math] changes in the range from [math]1[/math] to [math]i - 1[/math] and takes all integer values from this range.

Arguments of this operation:

  • [math]a[/math]:
    • for [math]p = 1[/math]: the element [math]a_{ji}[/math];
    • for [math]p \gt 1[/math]: the result of the operation corresponding to the vertex of the third group with coordinates [math]i, j, p - 1[/math];
  • [math]b[/math]: the result of the operation corresponding to the vertex of the second group with coordinates [math]p, i[/math];
  • [math]c[/math]: the result of the operation corresponding to the vertex of the second group with coordinates [math]p, j[/math].

The result of this operation is intermediate for the Cholesky algorithm.

The above graph is illustrated in Figs. 1 and 2 for [math]n = 4[/math]. In these figures, the vertices of the first group are highlighted in yellow and are marked by the letters SQ; the vertices of the second group are highlighted in green and are marked by the division sign; the vertices of the third group are highlighted in red and are marked by the letter F. The vertices corresponding to the results of operations (output data) are marked by large circles. The arcs doubling one another are depicted as a single one. The representation of the graph shown in Fig.1 corresponds to the classical definition. The graph of Fig.2 contains additional vertices corresponding to input data (blue) and to output data (pink).

Figure 1. The graph of the Cholesky algorithm without input and output data: SQ is the square-root operation, F is the operation a-bc, Div is division.
Figure 2. The graph of the Cholesky algorithm with input and output data: SQ is the square-root operation, F is the operation a-bc, Div is division, In and Out indicate input and output data.

1.8 Parallelization resource of the algorithm

In order to decompose a matrix of order [math]n[/math], a parallel version of the Cholesky algorithm should serially perform the following layers of operations:

  • [math]n[/math] layers for square roots (a single square root on each layer);
  • [math]n - 1[/math] layers for divisions (on each of the layers, the number of divisions is linear and, depending on a particular layer, varies from [math]1[/math] to [math]n - 1[/math]);
  • [math]n - 1[/math] layers of multiplications and [math]n - 1[/math] layers of addition/subtraction (on each of the layers, the number of these operations is quadratic and varies from [math]1[/math] to [math]\frac{n^2 - n}{2}[/math]).

Contrary to a serial version, in a parallel version the square-root and division operations require a significant part of overall computational time. The existence of isolated square roots on some layers of the parallel form may cause other difficulties for particular parallel computing architectures. In the case of programmable logic devices (PLD), for example, the operations of division, multiplication and addition/subtraction can be conveyorized, which is efficient in resource saving for programmable boards, whereas the isolated square roots acquire resources of such boards and force them to be out of action for a significant period of time. In the case of symmetric linear systems, the Cholesky decomposition is preferable compared to Gaussian elimination because of the reduction in computational time by a factor of two. However, this is not true in the case of its parallel version.

In addition, we should mention the fact that the accumulation mode requires multiplications and subtraction in double precision. In a parallel version, this means that almost all intermediate computations should be performed with data given in their double precision format. Contrary to a serial version, hence, this almost doubles the memory expenditure.

Thus, the Cholesky decomposition belongs to the class of algorithms of linear complexity in the sense of the height of its parallel form, whereas its complexity is quadratic in the sense of the width of its parallel form.

1.9 Input and output data of the algorithm

Input data: a dense matrix [math]A[/math] of order [math]n[/math]; its elements are denoted by [math]a_{ij}[/math]. The additional conditions:

  • [math]A[/math] is a symmetric matrix (i.e., [math]a_{ij}= a_{ji}, i, j = 1, \ldots, n[/math]).
  • [math]A[/math] is a positive definite matrix (i.e., [math]\vec{x}^T A \vec{x} \gt 0[/math] for any nonzero vector [math]\vec{x}[/math] of length [math]n[/math].

Amount of input data: [math]\frac{n (n + 1)}{2}[/math] ; since the matrix is symmetric, it is sufficient to store only its main diagonal and the elements above or under this diagonal. In practice, this storage saving scheme can be implemented in various ways. In the Numerical Analysis Library developed at the Moscow University Research Computing Center, for example, the lower triangle of the matrix [math]A[/math] is stored row-wise in a one-dimensional array of length [math]\frac{n (n + 1)}{2}[/math].

Output data: the lower triangular matrix [math]L[/math] whose elements are denoted by [math]l_{ij}[/math].

Amount of output data: [math]\frac{n (n + 1)}{2}[/math] ; since the matrix [math]L[/math] is lower triangular, it is sufficient to store only its main diagonal and the elements under this diagonal. In practice, this storage saving scheme can be implemented in various ways. In the above-mentioned library, for example, the matrix [math]L[/math] is stored row-wise in a one-dimensional array of length [math]\frac{n (n + 1)}{2}[/math].

1.10 Properties of the algorithm

In the case of unlimited computer resources, the ratio of the serial complexity to the parallel complexity is quadratic.

The computational power of the Cholesky algorithm considered as the ratio of the number of operations to the amount of input and output data is only linear.

The Cholesky is almost completely deterministic, which is ensured by the uniqueness theorem for this particular decomposition. Another order of associative operations may lead to the accumulation of round-off errors; however, the effect of this accumulation is not so large as in the case of not using the accumulation mode when computing dot products.

The information graph arcs from the vertices corresponding to the square-root and division operations can be considered as groups of data such that the function relating the multiplicity of these vertices and the number of these operations is a linear function of the matrix order and the vertex coordinates. These groups may contain «long» arcs; the remaining arcs are local.

The most known is the compact packing of a graph in the form of its projection onto the matrix triangle whose elements are recomputed by the packed operations. The «long» arcs can be eliminated and replaced by a combination of short-range arcs.

The following inequality is valid for the equivalent perturbation [math]M[/math] of the Cholesky decomposition: [math] ||M||_{E} \leq 2||\delta A||_{E}, [/math] where [math]\delta A[/math] is the perturbation of the matrix [math]A[/math] caused by the representation of the matrix elements in the computer memory. Such a slow growth of matrix elements during decomposition is due to the fact that the matrix is symmetric and positive definite.

2 Software implementation of the algorithm

2.1 Implementation peculiarities of the serial algorithm

In its simplest version without permuting the summation, the Cholesky decomposition can be represented in Fortran as 

	DO  I = 1, N
		S = A(I,I)
		DO  IP=1, I-1
			S = S - DPROD(A(I,IP), A(I,IP))
		END DO
		A(I,I) = SQRT (S)
		DO  J = I+1, N
			S = A(J,I)
			DO  IP=1, I-1
				S = S - DPROD(A(I,IP), A(J,IP))
			END DO
			A(J,I) = S/A(I,I)
		END DO
	END DO

Here [math]S[/math] is a double precision variable in the case of the accumulation mode.

A block version of the Cholesky algorithm is usually implemented in such a way that the scalar operations in its serial versions are replaced by the corresponding block-wise operations instead of using the loop unrolling and reordering techniques.

In order to ensure the locality of memory access in the Cholesky algorithm, in its Fortran implementation the original matrix and its decomposition are stored in the upper triangle instead of the lower triangle. The efficiency of such a version can be explained by the fact that Fortran stores matrices by columns and, hence, the computer programs in which the inner loops go up or down a column generate serial access to memory, contrary to the non-serial access when the inner loop goes across a row. This column orientation provides a significant improvement on computers with paging and cache memory.

There exists the following dot version of the Cholesky decomposition: the elements of the matrix [math]L[/math] are used as arguments in subsequent operations as soon as they are computed. This version can be illustrated as follows: 

	DO  I = 1, N
		A(I,I) = SQRT (A(I, I))
		DO  J = I+1, N
			A(J,I) = A(J,I)/A(I,I)
		END DO
		DO  K=I+1,N
			DO J = K, N
				A(J,K) = A(J,K) - A(J,I)*A(K,I)
			END DO	
		END DO
	END DO

As can be seen from the above program fragment, the array to store the original matrix and the output data should be declared as double precision for the accumulation mode. Note that the graph of the algorithm for this fragment and for the previous one is almost the same (the only distinction is that the DPROD function is used instead of multiplications).

2.2 Possible methods and considerations for parallel implementation of the algorithm

An analysis of the algorithm’s structure allows one to conclude that a large variety of its parallel implementations can be proposed. In the second version (see «Implementation peculiarities of the serial algorithm»), for example, all the inner loops are parallel, whereas, in the first version, only the inner loop over [math]J[/math] is parallel. Nevertheless, a simple parallelization technique causes a large number of data transfer between the processors at each step of the outer loop; this number is almost comparable with the number of arithmetic operations. Hence, it is reasonable to partition the computations into blocks with the corresponding partitioning of the data arrays before the allocation of operations and data between the processors of the computing system in use.

A good-enough decision depends on the features of a particular computer system. If the nodes of a multiprocessor computer are equipped with conveyors, it is reasonable to compute several dot products at once in parallel. Such a possibility exists in the case of programmable logic devices; in this case, however, the arithmetic speed is limited by a slow serial square-root operation. In principle, it is possible to apply the so-called «skew» parallelism; however, this approach is not used in practice because of complexity in the control structure of the algorithm’s implementation.

The dot version of the Cholesky algorithm is implemented in the basic libraries produced in Russia and in the libraries produced in western countries (for example, in LINPACK, LAPACK, ScaLAPACK, and others).

In the Russian libraries, as a rule, the accumulation mode is implemented to reduce the effect of round-off errors. In order to save computer memory, a packed representation of the matrices [math]A[/math] and [math]L[/math] is also used in the form of one-dimensional arrays.

In the modern western libraries, the dot Cholesky implementations are based on its LINPACK implementation utilizing the BLAS library. The dot product is computed in BLAS with no accumulation mode, but this mode can easily be implemented by minor modifications of the SDOT subprogram included in the BLAS library.

It is interesting to note that the original Cholesky decomposition is available in the largest western libraries, whereas the faster [math]LU[/math] decomposition algorithm without square-root operations is used only in special cases (for example, for tridiagonal matrices) when the number of diagonal elements is comparable with the number of off-diagonal elements.

2.3 Run results

2.4 Conclusions for different classes of computer architecture

Analyzing the performance of the ScaLAPACK library on supercomputers, we can conclude that, for large values of the matrix order [math]n[/math], the data exchanges reduce the execution efficiency, but to a smaller extent than the nonoptimality in the organization of computations within a single node. At the first stages, hence, it is necessary to optimize not a block algorithm but the subroutines used on individual processors, such as the dot Cholesky decomposition, matrix multiplications, etc. A number of possible directions of such an optimization are discussed below.

Generally speaking, the efficiency of the Cholesky algorithm cannot be high for parallel computer architectures. This fact can explained by the following property of its information structure: the square-root operations are the bottleneck of the Cholesky algorithm in comparison with the division operations or with the [math]a - bc[/math] operations, since these operations can easily be conveyorized or distributed among the nodes. In order to enhance the performance efficiency on parallel computers, hence, it is reasonable to use not the original Cholesky algorithm but its well-known modification without square-root operations; see the matrix decomposition in the form [math]L D L^T[/math][12].

3 References

  1. ↑ Commandant Benoit, Note sur une méthode de résolution des équations normales provenant de l'application de la méthode des moindres carrés à un système d'équations linéaires en nombre inférieur à celui des inconnues (Procédé du Commandant Cholesky), Bulletin Géodésique 2 (1924), 67-77.
  2. ↑ Banachiewicz T. Principles d'une nouvelle technique de la méthode des moindres carrês. Bull. Intern. Acad. Polon. Sci. A., 1938, 134-135.
  3. ↑ Banachiewicz T. Méthode de résoltution numérique des équations linéaires, du calcul des déterminants et des inverses et de réduction des formes quardatiques. Bull. Intern. Acad. Polon. Sci. A., 1938, 393-401.
  4. ↑ Voevodin V.V. Computational foundations of linear algebra Moscow: Nauka, 1977.
  5. ↑ Voevodin V.V., Kuznetsov Yu.A. Matrices and computations, Moscow: Nauka, 1984.
  6. ↑ Faddeev D.K., Faddeeva V.N. Computational methods in linear algebra. Moscow-Leningrad: Fizmatgiz, 1963.
  7. ↑ Kaporin I.E. High quality preconditioning of a general symmetric positive definite matrix based on its UTU + UTR + RTU-decomposition. Numer. Lin. Algebra Appl. (1998) Vol. 5, No. 6, 483-509.
  8. ↑ Kaporin I.E., Kon'shin I.N. Parallel solution of symmetric positive definite systems based on decomposition into overlapping blocks. Zhurn. Vychisl. Matem. i Matem. Fiziki. (2001) Vol. 41, No. 4, 515–528.
  9. ↑ Voevodin V.V. Mathematical foundations of parallel computing. Moscow: Moscow University Press, 1991. 345 pp.
  10. ↑ Voevodin V.V., Voevodin Vl.V. Parallel computing. St. Petersburg : BHV-Petersburg, St. 2002. 608 pp.
  11. ↑ Frolov A.V. Design principles and description of the Sigma language. Preprint N 236. Moscow: Institute of Numerical Mathematics, 1989.
  12. ↑ Krishnamoorthy A., Menon D. Matrix Inversion Using Cholesky Decomposition. 2013. eprint arXiv:1111.4144
Retrieved from "https://algowiki-project.org/algowiki/en/index.php?title=Cholesky_decomposition&oldid=10487"