It is often required to solve systems of linear equations with a symmetric positive-definite coefficient matrix. Of course, these systems could be solved by using algorithm on the basis of the LU decomposition as systems with a general matrix. However, taking into account the coefficient matrix properties we can use not the LU decomposition but Cholesky decomposition.

Both algorithms find the same solution. However the algorithm on the basis of the Cholesky decomposition is approximately two times as fast, so in case of a symmetric positive-definite coefficient matrix, it is reasonable to save time by using a specialized algorithm.

Note #1
Although the number of arithmetic operations in the Cholesky decomposition is twice as small as in the LU decomposition, the Cholesky decomposition program execution times may not be exactly two times as fast. The algorithms process matrices differently and depending on the way the matrix is stored in the main memory, algebraically equivalent operations could be executed in a different time subject to their carrying out over the rows or columns.

Note #2
Algorithm of a symmetric positive definite matrix condition number estimate can be used to estimate the matrix condition number.

It should be noted that the algorithm can solve systems with a symmetric positive-definite matrix only because the Cholesky decomposition exists only for such matrices. To solve systems of linear equations with an arbitrary symmetric matrix, use the algorithm on the basis of LDLT decomposition.

Subroutine description

As for the algorithm which is based on the LU decomposition, there are two subroutines for the algorithm. The first subroutine, SPDMatrixSolve, solves a system of linear equations with a symmetric positive-definite matrix which is given by its upper or lower triangle. This subroutine calls the SPDMatrixCholesky subroutine to perform the Cholesky decomposition, and then solves the obtained system by using the same method as when performing the LU decomposition.

The second subroutine, SPDMatrixCholeskySolve, is used to solve systems of linear equations whose Cholesky decomposition has already been performed. As with the LU decomposition, it is possible to solve a set of systems of linear equations with different right sides but with the same coefficient matrices because the Cholesky decomposition requires O(N ³) operations, and solving a system of linear equations with a coefficient matrix given by the Cholesky decomposition requires O(N ²) operations.

This algorithm is transferred from the LAPACK library.

Report a bug

C#

C# 1.0 source.
spdsolve.csharp.zip - Solution of the system with symmetric positive definite matrix

C++

C++ source.
spdsolve.cpp.zip - Solution of the system with symmetric positive definite matrix
ablas.zip - optimized basic linear algebra subroutines with SSE2 support (for C++ sources only)

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.
spdsolve.mpfr.zip - Solution of the system with symmetric positive definite matrix
mpfr.zip - precompiled Win32 MPFR/GMP binaries

Delphi

Delphi source.
Can be compiled under FPC (in Delphi compatibility mode).
spdsolve.delphi.zip - Solution of the system with symmetric positive definite matrix

Visual Basic 6

Visual Basic 6 source.
spdsolve.vb6.zip - Solution of the system with symmetric positive definite matrix

Zonnon ^beta

Zonnon source.
Zonnon is an experimental language developed at ETH Zurich.
See www.zonnon.ethz.ch for more information.
spdsolve.zonnon.zip - Solution of the system with symmetric positive definite matrix