Systems of Linear Equations



Choleski Decomposition

Choleski Decomposition

Choleski's method decomposes an n × n positive definite symmetric matrix A into the product of an n×n lower triangular matrix L and an n×n upper triangular matrix U, where
L = UT, i.e. A = L LT = UT U.  Choleski's algorithm proceeds as follows:
Evaluate the following pair of expressions for k = 0, . . . , n-1;
Lkk = [ Akk - (Lk0² + · · · + Lk,k-1² ) ]½,
and
Lik = [ Aik - (Li0 Lk0 + · · · + Li,k-1 Lk,k-1) ] / Lkk,   for i = k+1, . . . , n-1.
Note that Lik = 0 for k > i.  The matrix A = L LT.
After Choleski decomposition of A is performed, the solution to the system of linear equations A x = L LT x = B is solved by solving the system of linear equations L y = B by forward substitution for y, and then solving the system of linear equations LT x = y by backward substitution for x.

In case the actual inverse of A is needed, a function is provided to calculate A-1 = L-T L-1 using the Choleski decomposition of A.

Source code is provided for the six different versions of Choleski's L LT decomposition, the different versions reflect alternative ways of storing a positive definite symmetric matrix in computer memory.

An alternative to Choleski's method which avoids taking a square root decomposes an n×n positive definite symmetric matrix A into the product of an n×n unit lower triangular matrix L, an n×n diagonal matrix D, and an n×n unit upper triangular matrix U, where
L = UT, i.e. A = L D LT = UT D U.  Choleski's alternative algorithm proceeds as follows:
Evaluate the following pair of expressions for i = 0, . . . , n-1;
Lik dk = [ Aik - (Li0 d0 Lk0 + · · · + Li,k-1 dk-1 Lk,k-1) ]
and
di = [ Aii - (Li0 d0 Li0 + · · · + Lk,k-1 dk-1 Li,k-1 ) ],   for k = 0, . . . , i-1.
Note that Lik = 0 for k > i and Lkk = 1, for k = 0, . . . , n-1.  The matrix A = L D LT.
After Choleski decomposition of A is performed, the solution to the system of linear equations A x = L D LT x = B is solved by solving the system of linear equations L y = B by forward substitution for y, then solving D z = y for z by dividing yi by di, and then solving the system of linear equations LT x = z by backward substitution for x.

In case the actual inverse of A is needed, a function is provided to calculate A-1 = L-T D-1 L-1 using the Choleski decomposition of A.

Source code is provided for the six different versions of Choleski's L D LT decomposition, the different versions reflect alternative ways of storing a positive definite symmetric matrix in computer memory.

Function List

  • int Choleski_LU_Decomposition(double* A, int n)

    Given the postive definite symmetric n×n matrix A, the function Choleski_LU_Decomposition uses Choleski's algorithm to decompose A into the product of a lower triangular matrix and its transpose. The lower triangular part of the lower triangular matrix factor is returned in the lower triangular part of the A and the transpose of the lower triangular matrix factor is returned in the upper triangular part of the A so that the original matrix is replaced with the triangular factorizations.  Choleski_LU_Decomposition returns 0 if the decomposition was successful and returns -1 if the matrix is not positive definite symmetric to working accuracy.

  • int Choleski_LU_Solve(double *LU, double* B, double* x, int n)

    The function Choleski_LU_Solve solves the system of linear equations LU x = B for x, where LU is the LU decomposition of A evaluated by Choleski_LU_Decomposition.  Choleski_LU_Solve returns 0 if the solution was found and returns -1 if the matrix is singular.

  • int Choleski_LU_Inverse(double *LU, int n)

    The function Choleski_LU_Inverse computes the inverse of the A using the decomposition computed by Choleski_LU_Decomposition. The inverse of A = L U replaces the matrix argument LU. Choleski_LU_Inverse returns 0 if the inverse was computed and returns -1 if the matrix is singular.

  • int Choleski_LU_Decomposition_lt(double* A, int n)

    Given the positive definite symmetric n×n matrix A, due to the symmetry of A one need only store either the lower triangular part of A or the upper triangular part of A. For this routine it is assumed that only the lower triangular part of A is stored, i.e. A points to a memory location which contains A[0][0] immediately followed in memory by A[1][0], A[1][1], … ,A[n-1][0], A[n-1][1], … ,A[n-1][n-1]. Choleski_LU_Decomposition_lt uses Choleski's algorithm to decompose A into the product of a lower triangular matrix and its transpose. The lower triangular part of the lower triangular matrix factor is returned in the memory addresses previously occupied by A.  Choleski_LU_Decomposition_lt returns 0 if the decomposition was successful and returns -1 if the matrix is not positive definite symmetric to working accuracy.

  • int Choleski_LU_Solve_lt(double *LU, double* B, double* x, int n)

    Choleski_LU_Solve_lt solves the system of linear equations LU x = B for x, where LU is the LU decomposition of A returned from Choleski_LU_Decomposition_lt.  Choleski_LU_Solve_lt returns 0 if the solution was found and returns -1 if the matrix is singular.

  • int Choleski_LU_Inverse_lt(double *LU, int n)

    Choleski_LU_Inverse_lt computes the inverse of the A using the decomposition computed by Choleski_LU_Decomposition_lt. Only the lower triangular part of the inverse is returned in the memory locations previously occupied by the lower triangular part of the lower trangular matrix computed by Choleski_LU_Decomposition_lt. Choleski_LU_Inverse_lt returns 0 if the inverse was computed and returns -1 if the matrix is singular.

  • int Choleski_LU_Decomposition_ut(double* A, int n)

    Given the positive definite symmetric n×n matrix A, due to the symmetry of A one need only store either the lower triangular part of A or the upper triangular part of A. For this routine it is assumed that only the upper triangular part of A is stored, i.e. A points to a memory location which contains A[0][0] immediately followed in memory by A[0][1], … , A[0][n-1], A[1][1], … ,A[1][n-1], … ,A[n-1][n-1]. Choleski_LU_Decomposition_ut uses Choleski's algorithm to decompose A into the product of a lower triangular matrix and its transpose. The upper triangular part of the transpose of the lower triangular matrix factor is returned in the memory addresses previously occupied by A.  Choleski_LU_Decomposition_ut returns 0 if the decomposition was successful and returns -1 if the matrix is not positive definite symmetric to working accuracy.

  • int Choleski_LU_Solve_ut(double *LU, double* B, double* x, int n)

    Choleski_LU_Solve_ut solves the system of linear equations LU x = B for x, where LU is the LU decomposition of A returned from Choleski_LU_Decomposition_ut.  Choleski_LU_Solve_ut returns 0 if the solution was found and returns -1 if the matrix is singular.

  • int Choleski_LU_Inverse_ut(double *LU, int n)

    Choleski_LU_Inverse_ut computes the inverse of the A using the decomposition computed by Choleski_LU_Decomposition_ut. Only the upper triangular part of the inverse is returned in the memory locations previously occupied by the upper triangular part of the transpose of the lower trangular matrix computed by Choleski_LU_Decomposition_ut. Choleski_LU_Inverse_ut returns 0 if the inverse was computed and returns -1 if the matrix is singular.

  • int Choleski_LU_Decomposition_Band(double* A, int n, int bandwidth)

    A band matrix is a matrix whose elements vanish outsize a band along the main diagonal. For an arbitrary band matrix A, the bandwidth is defined as that number m such that if | i - j | > ( m - 1 ) / 2, then aij = 0.  For example a diagonal matrix has bandwidth of 1, a tridiagonal matrix has bandwidth of 3.

    Given the postive definite symmetric n×n band matrix A with bandwidth bandwidth, Choleski_LU_Decomposition_Band uses Choleski's algorithm to decompose A into the product of a lower triangular matrix and its transpose. The lower triangular part of the lower triangular matrix factor is returned in the lower triangular part of the A and the transpose of the lower triangular matrix factor is returned in the upper triangular part of the A so that the original matrix is replaced with the triangular factorizations.  Choleski_LU_Decomposition_Band returns 0 if the decomposition was successful and returns -1 if the matrix is not positive definite symmetric to working accuracy.

  • int Choleski_LU_Solve_Band(double *LU, double* B, double* x, int n, int bandwidth)

    Choleski_LU_Solve_Band solves the system of linear equations LU x = B for x, where LU is the LU decomposition of A returned from Choleski_LU_Decomposition_Band.  Choleski_LU_Solve_Band returns 0 if the solution was found and returns -1 if the matrix is singular.

    In general the inverse of a band matrix is not a band matrix, there is not algorithm presented here to calculate the inverse of a positive definite band matrix, however one can use Choleski_LU_Inverse.

  • int Choleski_LU_Decomposition_Band_band(double* A, int n, int bandwidth)

    A band matrix is a matrix whose elements vanish outsize a band along the main diagonal. For an arbitrary band matrix A, the bandwidth is defined as that number m such that if | i - j | > ( m - 1 ) / 2, then aij = 0.  For example a diagonal matrix has bandwidth of 1, a tridiagonal matrix has bandwidth of 3.

    For a band matrix only the elements of the band need to be stored in computer memory; therefore for each row, only the elements of the band are stored. The center of the band always contains the diagonal element of the matrix, elements of the band which do not correspond to a matrix element are set to 0.0. A matrix for which each row only the elements of the band are stored is said to be stored in band form.

    Given the postive definite symmetric n×n band matrix A with bandwidth bandwidth stored in band matrix form, Choleski_LU_Decomposition_Band_band uses Choleski's algorithm to decompose A into the product of a lower triangular matrix factor and its transpose. The original matrix is replaced with the triangular factorizations: the lower triangular part of the lower triangular matrix is returned in the lower triangular part of A and the transpose in the upper triangular part of A which taken together both factors are stored in band form.  Choleski_LU_Decomposition_Band_band returns 0 if the decomposition was successful and returns -1 if the matrix is not positive definite symmetric to working accuracy.

  • int Choleski_LU_Solve_Band_band(double *LU, double* B, double* x, int n, int bandwidth)

    Choleski_LU_Solve_Band_band solves the system of linear equations LU x = B for x, where LU is the LU decomposition of A returned from Choleski_LU_Decomposition_Band_band.  Choleski_LU_Solve_Band_band returns 0 if the solution was found and returns -1 if the matrix is singular.

    In general the inverse of a band matrix is not a band matrix, there is not algorithm presented here to calculate the inverse of a positive definite band matrix stored in band matrix form.

  • int Choleski_LU_Decomposition_Band_sband(double* A, int n, int bandwidth)

    A band matrix is a matrix whose elements vanish outsize a band along the main diagonal. For a symmetric band matrix A, the bandwidth is defined as that number m such that if | i - j | ≥ m, then aij = 0.  For example a diagonal matrix has bandwidth of 1, a symmetric tridiagonal matrix has bandwidth of 2.

    For a symmetric band matrix only the elements of the band to the left of the diagonal and the diagonal need to be stored in computer memory; therefore for each row, only the elements of the band to the left of the diagonal and the diagonal are stored. The right-most element of the band always contains the diagonal element of the matrix, elements of the band which do not correspond to a matrix element are set to 0.0. A symmetric matrix for which each row only the elements of the band to the left of the diagonal and the diagonal are stored is said to be stored in symmetric band form.

    Given the postive definite symmetric n×n band matrix A with bandwidth bandwidth which is stored in symmetric band matrix form Choleski_LU_Decomposition_Band_sband uses Choleski's algorithm to decompose A into the product of a lower triangular matrix and its transpose. The lower triangular part of the lower triangular matrix is returned in the lower triangular part of A and the transpose in the upper triangular part of A.  Choleski_LU_Decomposition_Band_sband returns 0 if the decomposition was successful and returns -1 if the matrix is not positive definite symmetric to working accuracy.

  • int Choleski_LU_Solve_Band_sband(double *LU, double* B, double* x, int n, int bandwidth)

    Choleski_LU_Solve_Band_sband solves the system of linear equations LU x = B for x, where LU is the LU decomposition of A returned from Choleski_LU_Decomposition_Band_sband.  Choleski_LU_Solve_Band_sband returns 0 if the solution was found and returns -1 if the matrix is singular.

    In general the inverse of a band matrix is not a band matrix, there is not algorithm presented here to calculate the inverse of a positive definite band matrix stored in band matrix form.

  • int Choleski_LDU_Decomposition(double* A, int n)

    Given the postive definite symmetric n × n matrix A, Choleski_LDU_Decomposition uses Choleski's algorithm to decompose A into the product of a unit lower triangular matrix, a diagonal matrix and the transpose of the unit lower triangular matrix. The non-diagonal part of the lower triangular part of the lower triangular matrix is returned in the non-diagonal lower triangular part of A, the diagonal of the diagonal matrix is stored in the diagonal part of A and the non-diagonal part of the transpose of the lower triangular matrix is stored in the non-diagonal part of the upper triangular part of A.  Choleski_LDU_Decomposition returns 0 if the decomposition was successful and returns -1 if the matrix is not positive definite symmetric to working accuracy.

  • int Choleski_LDU_Solve(double *LDU, double* B, double* x, int n)

    Choleski_LDU_Solve solves the system of linear equations LDU x = B for x, where LDU is the LDU decomposition of A returned from Choleski_LDU_Decomposition.  Choleski_LDU_Solve returns 0 if the solution was found and returns -1 if the matrix is singular.

  • int Choleski_LDU_Inverse(double *LDU, int n)

    Choleski_LDU_Inverse computes the inverse of the A using the decomposition computed by Choleski_LDU_Decomposition. Choleski_LDU_Inverse returns 0 if the inverse was computed and returns -1 if the matrix is singular.

  • int Choleski_LDU_Decomposition_lt(double* A, int n)

    Given the positive definite symmetric n × n matrix A, due to the symmetry of A we need only store either the lower triangular part of A or the upper triangular part of A. For this set of routines it is assumed that only the lower triangular part of A is stored, i.e. A points to a memory location which contains A[0][0] followed by A[1][0], A[1][1], … ,A[n-1][0], A[n-1][1], … ,A[n-1][n-1]. Choleski_LDU_Decomposition_lt uses Choleski's algorithm to decompose A into the product of a lower triangular matrix and its transpose. The lower triangular part of the lower triangular matrix is returned in the memory addresses previously occupied by A.  Choleski_LDU_Decomposition_lt returns 0 if the decomposition was successful and returns -1 if the matrix is not positive definite symmetric to working accuracy.

  • int Choleski_LDU_Solve_lt(double *LDU, double* B, double* x, int n)

    Choleski_LDU_Solve_lt solves the system of linear equations LU x = B for x, where LDU is the LDU decomposition of A returned from Choleski_LDU_Decomposition_lt.  Choleski_LDU_Solve_lt returns 0 if the solution was found and returns -1 if the matrix is singular.

  • int Choleski_LDU_Inverse_lt(double *LDU, int n)

    Choleski_LDU_Inverse_lt computes the inverse of the A using the decomposition computed by Choleski_LDU_Decomposition_lt. Only the lower triangular part of the inverse is returned in the memory locations previously occupied by the lower triangular part of the lower trangular matrix computed by Choleski_LDU_Decomposition_lt. Choleski_LDU_Inverse_lt returns 0 if the inverse was computed and returns -1 if the matrix is singular.
     
  • int Choleski_LDU_Decomposition_ut(double* A, int n)

    Given the positive definite symmetric n×n matrix A, due to the symmetry of A we need only store either the lower triangular part of A or the upper triangular part of A. For this set of routines it is assumed that only the upper triangular part of A is stored, i.e. A points to a memory location which contains A[0][0] followed by A[0][1], … , A[0][n-1], A[1][1], … ,A[1][n-1], … ,A[n-1][n-1]. Choleski_LDU_Decomposition_ut uses Choleski's algorithm to decompose A into the product of a lower triangular matrix and its transpose. The upper triangular part of the transpose of the lower triangular matrix is returned in the memory addresses previously occupied by A.  Choleski_LDU_Decomposition_ut returns 0 if the decomposition was successful and returns -1 if the matrix is not positive definite symmetric to working accuracy.

  • int Choleski_LDU_Solve_ut(double *LDU, double* B, double* x, int n)

    Choleski_LDU_Solve_ut solves the system of linear equations LDU x = B for x, where LDU is the LDU decomposition of A returned from Choleski_LDU_Decomposition_ut.  Choleski_LDU_Solve_ut returns 0 if the solution was found and returns -1 if the matrix is singular.

  • int Choleski_LDU_Inverse_ut(double *LDU, int n)

    Choleski_LDU_Inverse_ut computes the inverse of the A using the decomposition computed by Choleski_LDU_Decomposition_ut. Only the upper triangular part of the inverse is returned in the memory locations previously occupied by the upper triangular part of the transpose of th lower trangular matrix computed by Choleski_LDU_Decomposition_ut. Choleski_LDU_Inverse_ut returns 0 if the inverse was computed and returns -1 if the matrix is singular.

  • int Choleski_LDU_Decomposition_Band(double* A, int n, int bandwidth)

    A band matrix is a matrix whose elements vanish outsize a band along the main diagonal. For an arbitrary band matrix A, the bandwidth is defined as that number m such that if | i - j | > ( m - 1 ) / 2, then aij = 0.  For example a diagonal matrix has bandwidth of 1, a tridiagonal matrix has bandwidth of 3.

    Given the postive definite symmetric n×n band matrix A with bandwidth bandwidth, Choleski_LDU_Decomposition_Band uses Choleski's algorithm to decompose A into the product of a lower triangular matrix and its transpose. The lower triangular part of the lower triangular matrix is returned in the lower triangular part of A and the transpose in the upper triangular part of A.  Choleski_LDU_Decomposition_Band returns 0 if the decomposition was successful and returns -1 if the matrix is not positive definite symmetric to working accuracy.

  • int Choleski_LDU_Solve_Band(double *LDU, double* B, double* x, int n, int bandwidth)

    Choleski_LDU_Solve_Band solves the system of linear equations LDU x = B for x, where LDU is the LDU decomposition of A returned from Choleski_LDU_Decomposition_Band.  Choleski_LDU_Solve_Band returns 0 if the solution was found and returns -1 if the matrix is singular.

    In general the inverse of a band matrix is not a band matrix, there is not algorithm presented here to calculate the inverse of a positive definite band matrix, however one can use Choleski_LDU_Inverse.

  • int Choleski_LDU_Decomposition_Band_band(double* A, int n, int bandwidth)

    A band matrix is a matrix whose elements vanish outsize a band along the main diagonal. For an arbitrary band matrix A, the bandwidth is defined as that number m such that if | i - j | > ( m - 1 ) / 2, then aij = 0.  For example a diagonal matrix has bandwidth of 1, a tridiagonal matrix has bandwidth of 3.

    For a band matrix only the elements of the band need to be stored in computer memory; therefore for each row, only the elements of the band are stored. The center of the band always contains the diagonal element of the matrix, elements of the band which do not correspond to a matrix element are set to 0.0.

    Given the postive definite symmetric n×n band matrix A with bandwidth bandwidth which is stored in band matrix form Choleski_LDU_Decomposition_Band_band uses Choleski's algorithm to decompose A into the product of a lower triangular matrix and its transpose. The lower triangular part of the lower triangular matrix is returned in the lower triangular part of A and the transpose in the upper triangular part of A.  Choleski_LDU_Decomposition_Band_band returns 0 if the decomposition was successful and returns -1 if the matrix is not positive definite symmetric to working accuracy.

  • int Choleski_LDU_Solve_Band_band(double *LDU, double* B, double* x, int n, int bandwidth)

    Choleski_LDU_Solve_Band_band solves the system of linear equations LDU x = B for x, where LDU is the LDU decomposition of A returned from Choleski_LDU_Decomposition_Band_band.  Choleski_LDU_Solve_Band_band returns 0 if the solution was found and returns -1 if the matrix is singular.

    In general the inverse of a band matrix is not a band matrix, there is not algorithm presented here to calculate the inverse of a positive definite band matrix stored in band matrix form.

  • int Choleski_LDU_Decomposition_Band_sband(double* A, int n, int bandwidth)

    A band matrix is a matrix whose elements vanish outsize a band along the main diagonal. For a symmetric band matrix A, the bandwidth is defined as that number m such that if | i - j | ≥ m, then aij = 0.  For example a diagonal matrix has bandwidth of 1, a symmetric tridiagonal matrix has bandwidth of 2.

    For a symmetric band matrix only the elements of the band to the left of the diagonal and the diagonal need to be stored in computer memory; therefore for each row, only the elements of the band to the left of the diagonal and the diagonal are stored. The right-most element of the band always contains the diagonal element of the matrix, elements of the band which do not correspond to a matrix element are set to 0.0.

    Given the postive definite symmetric n × n band matrix A with bandwidth bandwidth which is stored in symmetric band matrix form Choleski_LDU_Decomposition_Band_sband uses Choleski's algorithm to decompose A into the product of a lower triangular matrix and its transpose. The lower triangular part of the lower triangular matrix is returned in the lower triangular part of A and the transpose in the upper triangular part of A.  Choleski_LDU_Decomposition_Band_sband returns 0 if the decomposition was successful and returns -1 if the matrix is not positive definite symmetric to working accuracy.

  • int Choleski_LDU_Solve_Band_sband(double *LDU, double* B, double* x, int n, int bandwidth)

    Choleski_LDU_Solve_Band_sband solves the system of linear equations LDU x = B for x, where LDU is the LDU decomposition of A returned from Choleski_LDU_Decomposition_Band_sband.  Choleski_LDU_Solve_Band_sband returns 0 if the solution was found and returns -1 if the matrix is singular.

    In general the inverse of a band matrix is not a band matrix, there is not algorithm presented here to calculate the inverse of a positive definite band matrix stored in band matrix form.

C Source Code

  • The file, choleski.c, contains the routines Choleski_LU_Decomposition( ), Choleski_LU_Solve( ), and Choleski_LU_Inverse( ).

  • The file, choleski_lt.c, contains the routines Choleski_LU_Decomposition_lt( ), Choleski_LU_Solve_lt( ), and Choleski_LU_Inverse_lt( ).

  • The file, choleski_ut.c, contains the routines Choleski_LU_Decomposition_ut( ), Choleski_LU_Solve_ut( ), and Choleski_LU_Inverse_ut( ).

  • The file, choleski_band.c, contains the routines Choleski_LU_Decomposition_Band( ) and Choleski_LU_Solve_Band( ).

  • The file, choleski_band_band.c, contains the routines Choleski_LU_Decomposition_Band_band( ) and Choleski_LU_Solve_Band_band( ).

  • The file, choleski_band_sband.c, contains the routines Choleski_LU_Decomposition_Band_sband( ) and
    Choleski_LU_Solve_Band_sband( ).

  • The file, choleski_ldu.c, contains the routines Choleski_LDU_Decomposition( ), Choleski_LDU_Solve( ), and Choleski_LDU_Inverse( ).

  • The file, choleski_ldu_lt.c, contains the routines Choleski_LDU_Decomposition_lt( ), Choleski_LDU_Solve_lt( ), and Choleski_LDU_Inverse_lt( ).

  • The file, choleski_ldu_ut.c, contains the routines Choleski_LDU_Decomposition_ut( ), Choleski_LDU_Solve_ut( ), and Choleski_LDU_Inverse_ut( ).

  • The file, choleski_ldu_band.c, contains the routines Choleski_LDU_Decomposition_Band( ) and Choleski_LDU_Solve_Band( ).

  • The file, choleski_ldu_band_band.c, contains the routines Choleski_LDU_Decomposition_Band_band( ) and Choleski_LDU_Solve_Band_band( ).

  • The file, choleski_ldu_band_sband.c, contains the routines Choleski_LDU_Decomposition_Band_sband( ) and Choleski_LDU_Solve_Band_sband( ).