14.8 Modified Cholesky Decomposition

The modified Cholesky decomposition is suitable for solving systems A x = b where A is a symmetric indefinite matrix. Such matrices arise in nonlinear optimization algorithms. The standard Cholesky decomposition requires a positive definite matrix and would fail in this case. Instead of resorting to a method like QR or SVD, which do not take into account the symmetry of the matrix, we can instead introduce a small perturbation to the matrix A to make it positive definite, and then use a Cholesky decomposition on the perturbed matrix. The resulting decomposition satisfies

P (A + E) P^T = L D L^T

where P is a permutation matrix, E is a diagonal perturbation matrix, L is unit lower triangular, and D is diagonal. If A is sufficiently positive definite, then the perturbation matrix E will be zero and this method is equivalent to the pivoted Cholesky algorithm. For indefinite matrices, the perturbation matrix E is computed to ensure that A + E is positive definite and well conditioned.

Function: int gsl_linalg_mcholesky_decomp (gsl_matrix * A, gsl_permutation * p, gsl_vector * E)

This function factors the symmetric, indefinite square matrix A into the Modified Cholesky decomposition P (A + E) P^T = L D L^T. On input, the values from the diagonal and lower-triangular part of the matrix A are used to construct the factorization. On output the diagonal of the input matrix A stores the diagonal elements of D, and the lower triangular portion of A contains the matrix L. Since L has ones on its diagonal these do not need to be explicitely stored. The upper triangular portion of A is unmodified. The permutation matrix P is stored in p on output. The diagonal perturbation matrix is stored in E on output. The parameter E may be set to NULL if it is not required.

Function: int gsl_linalg_mcholesky_solve (const gsl_matrix * LDLT, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x)

This function solves the perturbed system (A + E) x = b using the Cholesky decomposition of A + E held in the matrix LDLT and permutation p which must have been previously computed by gsl_linalg_mcholesky_decomp.

Function: int gsl_linalg_mcholesky_svx (const gsl_matrix * LDLT, const gsl_permutation * p, gsl_vector * x)

This function solves the perturbed system (A + E) x = b in-place using the Cholesky decomposition of A + E held in the matrix LDLT and permutation p which must have been previously computed by gsl_linalg_mcholesky_decomp. On input, x contains the right hand side vector b which is replaced by the solution vector on output.

Function: int gsl_linalg_mcholesky_rcond (const gsl_matrix * LDLT, const gsl_permutation * p, double * rcond, gsl_vector * work)

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