## Cholesky decomposition

*Definition*

Cholesky decomposition, along with LU decomposition, is one of the most popular triangular factorizations. Cholesky decomposition of a symmetric positive definite (SPD) matrix (or Hermitian positive definite, HPD) is quite similar to LU decomposition: *A* is represented as *A = LL*^{ H} (or, that essentially the same, as *A = U*^{ H}U). However there are differences too. First, there is no pivoting (but factorization is stable). Second, instead of two matrices (*L* and *U*) we have only one matrix multiplied by itself. Finally, Chelesky decomposition require two times less operation than LU decomposition of SPD/HPD matrix of the same size.

**Note #1**

Actual performance gain may be even higher than two because pivoting required by LU decomposition results in some performance penalty.

*Applications*

Main applications of Cholesky decomposition are linear systems and SPD/HPD matrix inversion.

*ALGLIB implementation*

ALGLIB package implements recursive version of Cholesky decomposition similar to the ATLAS one. Decomposition of NxN matrix is reduced to a sequence of approximately 0.5Nx0.5N decompositions. Recursion is stopped when matrix size becomes small enough to fit in the smallest L1 cache of the modern CPUs. At this point, non-recursive version is called. This approach has the following advantages:

- It makes effective use of CPU cache without explicitly specifying its size. With each iteration matrices become smaller until they fits completely in cache. Such algorithms are called 'cache oblivious'.
- Recursive algorithm makes extensive use of ALGLIB BLAS - mostly of the Level 3 BLAS, the most optimized part of ALGLIB package. The matrices being multiplied are usually square (which is more suited for optimization than a low rank updates generated by LAPACK routines).

These features allows us to achieve the maximum performance possible with our compiler/language. However, different compilers and different programming languages have different limitations. As for linear algebra algorithms, C++ version of ALGLIB is the most efficient implementation. Implementations in other languages have significantly lower performance.

*Subroutines*

Cholesky decomposition is calculated by **SPDMatrixCholesky** and **HPDMatrixCholesky** subroutines. Because *A* is symmetric/Hermitian, only upper or lower triangle may be specified. Correspondingly, as result we will get either lower triangular *L* or upper triangular *U*.

*Manual entries*

*This article is intended for personal use only.*

## Download ALGLIB