Next: Tridiagonal Decomposition of Real Symmetric Matrices, Previous: Pivoted Cholesky Decomposition, Up: Linear Algebra [Index]

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`.