14.7 Pivoted Cholesky Decomposition

A symmetric, positive definite square matrix A has an alternate Cholesky decomposition into a product of a lower unit triangular matrix L, a diagonal matrix D and L^T, given by L D L^T. This is equivalent to the Cholesky formulation discussed above, with the standard Cholesky lower triangular factor given by L D^{1 \over 2}. For ill-conditioned matrices, it can help to use a pivoting strategy to prevent the entries of D and L from growing too large, and also ensure D_1 \ge D_2 \ge \cdots \ge D_n > 0, where D_i are the diagonal entries of D. The final decomposition is given by

P A P^T = L D L^T

where P is a permutation matrix.

Function: int gsl_linalg_pcholesky_decomp (gsl_matrix * A, gsl_permutation * p)

This function factors the symmetric, positive-definite square matrix A into the Pivoted Cholesky decomposition P A 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.

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

This function solves the system A x = b using the Pivoted Cholesky decomposition of A held in the matrix LDLT and permutation p which must have been previously computed by gsl_linalg_pcholesky_decomp.

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

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

Function: int gsl_linalg_pcholesky_decomp2 (gsl_matrix * A, gsl_permutation * p, gsl_vector * S)

This function computes the pivoted Cholesky factorization of the matrix S A S, where the input matrix A is symmetric and positive definite, and the diagonal scaling matrix S is computed to reduce the condition number of A as much as possible. See Cholesky Decomposition for more information on the matrix S. The Pivoted Cholesky decomposition satisfies P S A S 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 scaling transformation is stored in S on output.

Function: int gsl_linalg_pcholesky_solve2 (const gsl_matrix * LDLT, const gsl_permutation * p, 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 Pivoted Cholesky decomposition of S A S held in the matrix LDLT, permutation p, and vector S, which must have been previously computed by gsl_linalg_pcholesky_decomp2.

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

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

Function: int gsl_linalg_pcholesky_invert (const gsl_matrix * LDLT, const gsl_permutation * p, gsl_matrix * Ainv)

This function computes the inverse of the matrix A, using the Pivoted Cholesky decomposition stored in LDLT and p. On output, the matrix Ainv contains A^{-1}.

Function: int gsl_linalg_pcholesky_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 symmetric positive definite matrix A, using its pivoted Cholesky decomposition provided in LDLT. 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.