Contents
Main
Site map
Links
Forum
Site and author
News
Contact

Cholesky decomposition

Contents

    Definition
    Applications
    ALGLIB implementation
    Subroutines
    Manual entries

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 HU). 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

C++ trfac.h   
C# trfac.cs   
MPFR trfac.h   
Delphi trfac.pas   
FreePascal trfac.pas   
VBA trfac.bas   

This article is intended for personal use only.

Download ALGLIB

C#

C# source.

alglib-2.5.0.csharp.zip

 

C++

C++ source.

alglib-2.5.0.cpp.zip

 

C++, multiple precision arithmetic

C++ source. MPFR/GMP is used.

GMP source is available from gmplib.org. MPFR source is available from www.mpfr.org.

alglib-2.5.0.mpfr.zip

 

FreePascal

FreePascal source.

alglib-2.5.0.freepascal.zip

 

Delphi

Delphi source.

alglib-2.5.0.delphi.zip

 

Visual Basic

VBA source.

alglib-2.5.0.vb6.zip

 


 
 
Sergey Bochkanov, Vladimir Bystritsky
Copyright © 1999-2010