14.6 Cholesky Decomposition

A symmetric, positive definite square matrix A has a Cholesky decomposition into a product of a lower triangular matrix L and its transpose L^T,

A = L L^T

This is sometimes referred to as taking the square-root of a matrix. The Cholesky decomposition can only be carried out when all the eigenvalues of the matrix are positive. This decomposition can be used to convert the linear system A x = b into a pair of triangular systems (L y = b, L^T x = y), which can be solved by forward and back-substitution.

If the matrix A is near singular, it is sometimes possible to reduce the condition number and recover a more accurate solution vector x by scaling as

( S A S ) ( S^(-1) x ) = S b

where S is a diagonal matrix whose elements are given by S_{ii} = 1/\sqrt{A_{ii}}. This scaling is also known as Jacobi preconditioning. There are routines below to solve both the scaled and unscaled systems.

Function: int gsl_linalg_cholesky_decomp1 (gsl_matrix * A)

Function: int gsl_linalg_complex_cholesky_decomp (gsl_matrix_complex * A)

These functions factorize the symmetric, positive-definite square matrix A into the Cholesky decomposition A = L L^T (or A = L L^H for the complex case). On input, the values from the diagonal and lower-triangular part of the matrix A are used (the upper triangular part is ignored). On output the diagonal and lower triangular part of the input matrix A contain the matrix L, while the upper triangular part is unmodified. If the matrix is not positive-definite then the decomposition will fail, returning the error code GSL_EDOM.

When testing whether a matrix is positive-definite, disable the error handler first to avoid triggering an error.

Function: int gsl_linalg_cholesky_decomp (gsl_matrix * A)

This function is now deprecated and is provided only for backward compatibility.

Function: int gsl_linalg_cholesky_solve (const gsl_matrix * cholesky, const gsl_vector * b, gsl_vector * x)

Function: int gsl_linalg_complex_cholesky_solve (const gsl_matrix_complex * cholesky, const gsl_vector_complex * b, gsl_vector_complex * x)

These functions solve the system A x = b using the Cholesky decomposition of A held in the matrix cholesky which must have been previously computed by gsl_linalg_cholesky_decomp or gsl_linalg_complex_cholesky_decomp.

Function: int gsl_linalg_cholesky_svx (const gsl_matrix * cholesky, gsl_vector * x)

Function: int gsl_linalg_complex_cholesky_svx (const gsl_matrix_complex * cholesky, gsl_vector_complex * x)

These functions solve the system A x = b in-place using the Cholesky decomposition of A held in the matrix cholesky which must have been previously computed by gsl_linalg_cholesky_decomp or gsl_linalg_complex_cholesky_decomp. On input x should contain the right-hand side b, which is replaced by the solution on output.

Function: int gsl_linalg_cholesky_invert (gsl_matrix * cholesky)

Function: int gsl_linalg_complex_cholesky_invert (gsl_matrix_complex * cholesky)

These functions compute the inverse of a matrix from its Cholesky decomposition cholesky, which must have been previously computed by gsl_linalg_cholesky_decomp or gsl_linalg_complex_cholesky_decomp. On output, the inverse is stored in-place in cholesky.

Function: int gsl_linalg_cholesky_decomp2 (gsl_matrix * A, gsl_vector * S)

This function calculates a diagonal scaling transformation S for the symmetric, positive-definite square matrix A, and then computes the Cholesky decomposition S A S = L L^T. On input, the values from the diagonal and lower-triangular part of the matrix A are used (the upper triangular part is ignored). On output the diagonal and lower triangular part of the input matrix A contain the matrix L, while the upper triangular part of the input matrix is overwritten with L^T (the diagonal terms being identical for both L and L^T). If the matrix is not positive-definite then the decomposition will fail, returning the error code GSL_EDOM. The diagonal scale factors are stored in S on output.

When testing whether a matrix is positive-definite, disable the error handler first to avoid triggering an error.

Function: int gsl_linalg_cholesky_solve2 (const gsl_matrix * cholesky, const gsl_vector * S, const gsl_vector * b, gsl_vector * x)

This function solves the system (S A S) (S^{-1} x) = S b using the Cholesky decomposition of S A S held in the matrix cholesky which must have been previously computed by gsl_linalg_cholesky_decomp2.

Function: int gsl_linalg_cholesky_svx2 (const gsl_matrix * cholesky, const gsl_vector * S, gsl_vector * x)

This function solves the system (S A S) (S^{-1} x) = S b in-place using the Cholesky decomposition of S A S held in the matrix cholesky which must have been previously computed by gsl_linalg_cholesky_decomp2. On input x should contain the right-hand side b, which is replaced by the solution on output.

Function: int gsl_linalg_cholesky_scale (const gsl_matrix * A, gsl_vector * S)

This function calculates a diagonal scaling transformation of the symmetric, positive definite matrix A, such that S A S has a condition number within a factor of N of the matrix of smallest possible condition number over all possible diagonal scalings. On output, S contains the scale factors, given by S_i = 1/\sqrt{A_{ii}}. For any A_{ii} \le 0, the corresponding scale factor S_i is set to 1.

Function: int gsl_linalg_cholesky_scale_apply (gsl_matrix * A, const gsl_vector * S)

This function applies the scaling transformation S to the matrix A. On output, A is replaced by S A S.

Function: int gsl_linalg_cholesky_rcond (const gsl_matrix * cholesky, double * rcond, gsl_vector * work)

This function estimates the reciprocal condition number (using the 1-norm) of the symmetric positive definite matrix A, using its Cholesky decomposition provided in cholesky. The reciprocal condition number estimate, defined as 1 / (||A||_1 \cdot ||A^{-1}||_1), is stored in rcond. Additional workspace of size 3 N is required in work.