Cholesky decomposition

The Cholesky decomposition (also Cholesky factorization ) (after André-Louis Cholesky , 1875-1918) describes in numerical mathematics a decomposition of a symmetrical positive definite matrix into a product of a lower triangular matrix and its transpose . It was developed by Cholesky before 1914 as part of the triangulation of Crete by the Service géographique de l'armée . The concept can also be defined more generally for Hermitian matrices .

Areas of application

When using the least squares method , one possibility is to solve the minimization problems that arise via the normal equations, which have a symmetrical positive definite system matrix. This is possible with the help of the Cholesky decomposition and this was the motivation for Cholesky to develop the decomposition. With the Gauss-Newton method , a system of equations has to be solved for each iteration step, which can be determined with the Cholesky method.

The Cholesky decomposition can also be used to obtain a preconditioning method for linear systems of equations with a positive definite matrix; for this purpose there are especially the variants of the incomplete Cholesky decomposition and the modified incomplete Cholesky decomposition .

At the same time, the decomposition represents a test whether a given symmetrical matrix is ​​positive definite. Otherwise one of the entries on the main diagonal is negative, so that the root cannot be drawn, or equal , so that the entry cannot be divided. The algorithm stops in both cases. The Cholesky decomposition can also be used to determine the determinant of the matrix  because it holds . ${\displaystyle 0}$ ${\displaystyle A}$ ${\displaystyle \det A=\prod _{i=1}^{n}G_{ii}^{2}}$

Outside of mathematics, the Cholesky decomposition is also used in econometric research into macroeconomic relationships. In so-called vector autoregressive models (VAR), the sequence in which the endogenous variables are influenced is determined.

In addition, it is also used in the Monte Carlo simulation in order to bring predefined correlations into independently generated random number sequences (as a discretization of stochastic processes).

Formulation and application

Every symmetrical , positively definite matrix can be unique in the form ${\displaystyle A\in \mathbb {R} ^{n\times n}}$

${\displaystyle A=LDL^{T}}$

to be written. There is a standardized lower triangular matrix and a diagonal matrix with positive entries. With the square root of and the matrix factor  defined by ${\displaystyle L}$ ${\displaystyle D}$ ${\displaystyle D}$ ${\displaystyle G}$

${\displaystyle D=D^{1/2}D^{1/2}}$

and

${\displaystyle G:=LD^{1/2}}$,

the Cholesky decomposition - equivalent - is also formulated as

${\displaystyle A=GG^{T}}$.

If the Cholesky decomposition has been calculated, the system of equations can be solved efficiently by inserting it forwards and backwards : ${\displaystyle Ax=b}$

• By substituting forwards: Solving the linear system of equations ${\displaystyle Gy=b}$

• By then inserting it backwards: Solving the linear equation system ${\displaystyle G^{T}x=y.}$

calculation

If one sets , one obtains for the elements of : ${\displaystyle A=GG^{T}\in \mathbb {R} ^{n\times n}}$ ${\displaystyle A=\left(a_{ij}\right)_{ij}}$

${\displaystyle a_{ij}=\sum \limits _{k=1}^{j}g_{ik}g_{jk}\,,\quad i\geq j.}$

This relationship leads directly to the following formulas for : ${\displaystyle G=\left(g_{ij}\right)_{ij}}$

${\displaystyle {\begin{aligned}g_{ij}&={\begin{cases}0&\mathrm {f{\ddot {u}}r} \ i<j\\{\sqrt {a_{ii}-\sum \limits _{k=1}^{i-1}g_{ik}^{2}}}&\mathrm {f{\ddot {u}}r} \ i=j\\{\frac {1}{g_{jj}}}\left(a_{ij}-\sum \limits _{k=1}^{j-1}g_{ik}g_{jk}\right)&\mathrm {f{\ddot {u}}r} \ i>j.\end{cases}}\\g_{i1}&={\begin{cases}{\sqrt {a_{11}}}&\mathrm {f{\ddot {u}}r} \ i=1\\{\frac {a_{i1}}{g_{11}}}={\frac {a_{i1}}{\sqrt {a_{11}}}}&\mathrm {f{\ddot {u}}r} \ i>1.\end{cases}}\end{aligned}}}$

With this algorithm it is important to carry out the order of the calculated matrix entries correctly. The entries are calculated column by column and starting with the lowest row index.

The decomposition is calculated in the same way for and : ${\displaystyle A=LDL^{T}}$ ${\displaystyle L=\left(l_{ij}\right)_{ij}}$ ${\displaystyle D=\left(d_{ij}\right)_{ij}}$

${\displaystyle {\begin{aligned}d_{ij}&={\begin{cases}0&\mathrm {f{\ddot {u}}r} \ i\neq j\\a_{ii}-\sum _{k=1}^{i-1}l_{ik}^{2}d_{kk}&\mathrm {f{\ddot {u}}r} \ i=j\end{cases}}\\l_{ij}&={\begin{cases}0&\mathrm {f{\ddot {u}}r} \ i<j\\1&\mathrm {f{\ddot {u}}r} \ i=j\\{\frac {1}{d_{jj}}}\left(a_{ij}-\sum _{k=1}^{j-1}l_{ik}l_{jk}d_{kk}\right)&\mathrm {f{\ddot {u}}r} \;i>j.\end{cases}}\end{aligned}}}$

With these algorithms, too, it is important to choose the correct order of the calculated matrix entries. First you have to calculate the entry for the index and then the -th column of the matrix , i.e.: for . ${\displaystyle j=1,\dotsc ,n}$ ${\displaystyle d_{jj}}$ ${\displaystyle j}$ ${\displaystyle L}$ ${\displaystyle l_{ij}}$ ${\displaystyle i=j+1,\dotsc ,n}$

Effort and stability

The Cholesky decomposition is numerically stable . In comparison, the elimination method according to Gauss with its algorithmic implementation, the LR decomposition , requires about twice as many operations, since not only one matrix  , but two factors  and have to be calculated. The Cholesky decomposition  involves arithmetic operations, including  multiplications,  divisions and root operations. ${\displaystyle G}$ ${\displaystyle L}$ ${\displaystyle R}$ ${\displaystyle {\tfrac {1}{3}}n^{3}+O(n^{2})}$ ${\displaystyle {\tfrac {1}{6}}n^{3}}$ ${\displaystyle {\tfrac {1}{2}}n^{2}}$ ${\displaystyle n}$

Pseudocode

The calculations in the above formulas can be performed in a number of ways. The variant named after Tadeusz Banachiewicz calculates the lower triangular matrix line by line. In pseudocode , the procedure for breaking down the matrix into the form looks like this : ${\displaystyle A}$ ${\displaystyle GG^{T}}$

Accesses (white) and writes (yellow).

    For i = 1 To n
        For j = 1 To i
            Summe = a(i, j)
            For k = 1 To j-1
                Summe = Summe - a(i, k) * a(j, k)
            If i > j Then
                a(i, j) = Summe / a(j, j)   // Untere Dreiecksmatrix
            Else If Summe > 0 Then          // Diagonalelement
                a(i, i) = Sqrt(Summe)       // ... ist immer groesser Null
            Else
                ERROR                       // Die Matrix ist (wenigstens numerisch) nicht symmetrisch positiv definit

The running indices in the code correspond to the mathematical notation of elements of the matrix . Here, the number of rows and at the same time the number of columns of the matrix , are auxiliary variables and sum . The algorithm works in-place ; that is, it modifies the matrix to include the lower triangular matrix . There is therefore no need to allocate new storage space. ${\displaystyle i,j=1,\,,\ldots \,,n}$ ${\displaystyle A=\left(a_{ij}\right)_{ij}}$ ${\displaystyle n}$ ${\displaystyle A}$ ${\displaystyle k}$ ${\displaystyle A}$ ${\displaystyle G}$

The algorithm only processes the lower left triangular matrix of , the elements for do not need to be assigned values ​​(since the matrix is symmetrical according to the assumption), and if they contain values, these are not changed. If one searches for the Cholesky decomposition according to , then the matrix elements from above the diagonal ( ) are to be set equal to 0. ${\displaystyle A=\left(a_{ij}\right)_{ij}}$ ${\displaystyle a_{ij}}$ ${\displaystyle i<j}$ ${\displaystyle A}$ ${\displaystyle G}$ ${\displaystyle A=GG^{T}}$ ${\displaystyle A}$ ${\displaystyle i<j}$

Web links

• taramath online tool for calculating the Cholesky decomposition of symmetrical and positive definite matrices.

literature

• Hans Rudolf Schwarz, Norbert Köckler: Numerical Mathematics. 5th edition. Teubner, Stuttgart 2004, ISBN 3-519-42960-8 .
• Gene H. Golub, Charles F. Van Loan: Matrix computations. 3rd edition. Johns Hopkins University Press, 1996, ISBN 0-8018-5414-8 .
• Michael Saunders: Commentary - Major Cholesky Would Feel Proud. In: ORSA Journal on Computing , 6, 1994, pp. 23-27.

Individual evidence

1. ↑ Andreas Meister: Numerics of linear systems of equations. 5th edition. Vieweg, Wiesbaden 2015, ISBN 3-528-13135-7 , p. 49.