Next: Pivoted Cholesky Decomposition, Previous: Singular Value Decomposition, Up: Linear Algebra [Index]

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