## File: Modified-Cholesky-Decomposition.html

package info (click to toggle)
gsl-ref-html 2.3-1
• area: non-free
• in suites: bullseye, buster, sid
• size: 6,876 kB
• ctags: 4,574
• sloc: makefile: 35
 file content (146 lines) | stat: -rw-r--r-- 8,799 bytes parent folder | download
 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146  GNU Scientific Library – Reference Manual: Modified Cholesky Decomposition

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.