Cholesky factorization

It is well-known that Hermitian positive-definite (HPD) matrices can be decomposed into the form \(A = L L^H\) or \(A = U^H U\), where \(L=U^H\) is lower triangular, and Cholesky factorization provides such an \(L\) (or \(U\)) given an HPD \(A\). If \(A\) is found to be numerically indefinite, then a NonHPDMatrixException will be thrown.

Implementation

Subroutines

Test driver

Without pivoting

The following routines overwrite the uplo triangle of the Hermitian positive-definite matrix A with its Cholesky factor without performing any pivoting.

Factorization

C++ API

void Cholesky(UpperOrLower uplo, Matrix<F> &A)

void Cholesky(UpperOrLower uplo, AbstractDistMatrix<F> &A)

C API

ElError ElCholesky_s(ElUpperOrLower uplo, ElMatrix_s A)

ElError ElCholesky_d(ElUpperOrLower uplo, ElMatrix_d A)

ElError ElCholesky_c(ElUpperOrLower uplo, ElMatrix_c A)

ElError ElCholesky_z(ElUpperOrLower uplo, ElMatrix_z A)

ElError ElCholeskyDist_s(ElUpperOrLower uplo, ElDistMatrix_s A)

ElError ElCholeskyDist_d(ElUpperOrLower uplo, ElDistMatrix_d A)

ElError ElCholeskyDist_c(ElUpperOrLower uplo, ElDistMatrix_c A)

ElError ElCholeskyDist_z(ElUpperOrLower uplo, ElDistMatrix_z A)

Solving linear systems with the factorization

C++ API

void cholesky::SolveAfter(UpperOrLower uplo, Orientation orientation, const Matrix<F> &A, Matrix<F> &B)

void cholesky::SolveAfter(UpperOrLower uplo, Orientation orientation, const AbstractDistMatrix<F> &A, AbstractDistMatrix<F> &B)

C API

ElError ElSolveAfterCholesky_s(ElUpperOrLower uplo, ElOrientation orientation, ElConstMatrix_s A, ElMatrix_s B)

ElError ElSolveAfterCholesky_d(ElUpperOrLower uplo, ElOrientation orientation, ElConstMatrix_d A, ElMatrix_d B)

ElError ElSolveAfterCholesky_c(ElUpperOrLower uplo, ElOrientation orientation, ElConstMatrix_c A, ElMatrix_c B)

ElError ElSolveAfterCholesky_z(ElUpperOrLower uplo, ElOrientation orientation, ElConstMatrix_z A, ElMatrix_z B)

ElError ElSolveAfterCholeskyDist_s(ElUpperOrLower uplo, ElOrientation orientation, ElConstDistMatrix_s A, ElDistMatrix_s B)

ElError ElSolveAfterCholeskyDist_d(ElUpperOrLower uplo, ElOrientation orientation, ElConstDistMatrix_d A, ElDistMatrix_d B)

ElError ElSolveAfterCholeskyDist_c(ElUpperOrLower uplo, ElOrientation orientation, ElConstDistMatrix_c A, ElDistMatrix_c B)

ElError ElSolveAfterCholeskyDist_z(ElUpperOrLower uplo, ElOrientation orientation, ElConstDistMatrix_z A, ElDistMatrix_z B)

Full pivoting

The following routines perform Cholesky factorization with diagonal pivoting, which can be shown to be equivalent to full pivoting for Hermitian positive-definite matrices. On exit, the vector \(p\) consists of the nonzero column indices of the permutation matrix \(P\) such that either \(P A P^T = L L^H\) or \(P A P^T = U^H U\).

Factorization

C++ API

void Cholesky(UpperOrLower uplo, Matrix<F> &A, Matrix<int> &p)

void Cholesky(UpperOrLower uplo, AbstractDistMatrix<F> &A, AbstractDistMatrix<int> &p)

C API

ElError ElCholeskyPiv_s(ElUpperOrLower uplo, ElMatrix_s A, ElMatrix_i p)

ElError ElCholeskyPiv_d(ElUpperOrLower uplo, ElMatrix_d A, ElMatrix_i p)

ElError ElCholeskyPiv_c(ElUpperOrLower uplo, ElMatrix_c A, ElMatrix_i p)

ElError ElCholeskyPiv_z(ElUpperOrLower uplo, ElMatrix_z A, ElMatrix_i p)

ElError ElCholeskyPivDist_s(ElUpperOrLower uplo, ElDistMatrix_s A, ElMatrix_i p)

ElError ElCholeskyPivDist_d(ElUpperOrLower uplo, ElDistMatrix_d A, ElMatrix_i p)

ElError ElCholeskyPivDist_c(ElUpperOrLower uplo, ElDistMatrix_c A, ElMatrix_i p)

ElError ElCholeskyPivDist_z(ElUpperOrLower uplo, ElDistMatrix_z A, ElMatrix_i p)

Solving linear systems with the factorization

C++ API

void cholesky::SolveAfter(UpperOrLower uplo, Orientation orientation, const Matrix<F> &A, Matrix<F> &B, Matrix<int> &p)

void cholesky::SolveAfter(UpperOrLower uplo, Orientation orientation, const AbstractDistMatrix<F> &A, AbstractDistMatrix<F> &B, AbstractDistMatrix<int> &p)

C API

ElError ElSolveAfterCholeskyPiv_s(ElUpperOrLower uplo, ElOrientation orientation, ElConstMatrix_s A, ElConstMatrix_i p, ElMatrix_s B)

ElError ElSolveAfterCholeskyPiv_d(ElUpperOrLower uplo, ElOrientation orientation, ElConstMatrix_d A, ElConstMatrix_i p, ElMatrix_d B)

ElError ElSolveAfterCholeskyPiv_c(ElUpperOrLower uplo, ElOrientation orientation, ElConstMatrix_c A, ElConstMatrix_i p, ElMatrix_c B)

ElError ElSolveAfterCholeskyPiv_z(ElUpperOrLower uplo, ElOrientation orientation, ElConstMatrix_z A, ElConstMatrix_i p, ElMatrix_z B)

ElError ElSolveAfterCholeskyPivDist_s(ElUpperOrLower uplo, ElOrientation orientation, ElConstDistMatrix_s A, ElConstDistMatrix_i p, ElDistMatrix_s B)

ElError ElSolveAfterCholeskyPivDist_d(ElUpperOrLower uplo, ElOrientation orientation, ElConstDistMatrix_d A, ElConstDistMatrix_i p, ElDistMatrix_d B)

ElError ElSolveAfterCholeskyPivDist_c(ElUpperOrLower uplo, ElOrientation orientation, ElConstDistMatrix_c A, ElConstDistMatrix_i p, ElDistMatrix_c B)

ElError ElSolveAfterCholeskyPivDist_z(ElUpperOrLower uplo, ElOrientation orientation, ElConstDistMatrix_z A, ElConstDistMatrix_i p, ElDistMatrix_z B)

Low-rank updates to a factorization

The following routines update an existing Cholesky factorization to incorporate the modification \(\alpha V V^H\) to the original matrix. The current algorithm uses Householder transformations for updates (\(\alpha \ge 0\)) and hyperbolic Householder transformations for downdates.

C++ API

void CholeskyMod(UpperOrLower uplo, Matrix<F> &T, Base<F> &alpha, Matrix<F> &V)

void CholeskyMod(UpperOrLower uplo, AbstractDistMatrix<F> &T, Base<F> &alpha, AbstractDistMatrix<F> &V)

C API

ElError ElCholeskyMod_s(ElUpperOrLower uplo, ElMatrix_s T, float alpha, ElMatrix_s V)

ElError ElCholeskyMod_d(ElUpperOrLower uplo, ElMatrix_d T, double alpha, ElMatrix_d V)

ElError ElCholeskyMod_c(ElUpperOrLower uplo, ElMatrix_c T, float alpha, ElMatrix_c V)

ElError ElCholeskyMod_z(ElUpperOrLower uplo, ElMatrix_z T, double alpha, ElMatrix_z V)

ElError ElCholeskyModDist_s(ElUpperOrLower uplo, ElDistMatrix_s T, float alpha, ElDistMatrix_s V)

ElError ElCholeskyModDist_d(ElUpperOrLower uplo, ElDistMatrix_d T, double alpha, ElDistMatrix_d V)

ElError ElCholeskyModDist_c(ElUpperOrLower uplo, ElDistMatrix_c T, float alpha, ElDistMatrix_c V)

ElError ElCholeskyModDist_z(ElUpperOrLower uplo, ElDistMatrix_z T, double alpha, ElDistMatrix_z V)

Rank-deficient factorization

It is possible to compute the Cholesky factor of a Hermitian positive semi-definite (HPSD) matrix through its eigenvalue decomposition, though it is significantly more expensive than the HPD case: Let \(A = U \Lambda U^H\) be the eigenvalue decomposition of \(A\), where all entries of \(\Lambda\) are non-negative. Then \(B = U \sqrt \Lambda U^H\) is the matrix square root of \(A\), i.e., \(B B = A\), and it follows that the QR and LQ factorizations of \(B\) yield Cholesky factors of \(A\):

\[A = B B = B^H B = (Q R)^H (Q R) = R^H Q^H Q R = R^H R,\]

and

\[A = B B = B B^H = (L Q) (L Q)^H = L Q Q^H L^H = L L^H.\]

If \(A\) is found to have eigenvalues less than \(-n \epsilon \| A \|_2\), then a NonHPSDMatrixException will be thrown.

Example driver

The following routines overwrite the uplo triangle of the potentially singular matrix A with its Cholesky factor.

C++ API

void HPSDCholesky(UpperOrLower uplo, Matrix<F> &A)

void HPSDCholesky(UpperOrLower uplo, AbstractDistMatrix<F> &A)

C API

ElError ElHPSDCholesky_s(ElUpperOrLower uplo, ElMatrix_s A)

ElError ElHPSDCholesky_d(ElUpperOrLower uplo, ElMatrix_d A)

ElError ElHPSDCholesky_c(ElUpperOrLower uplo, ElMatrix_c A)

ElError ElHPSDCholesky_z(ElUpperOrLower uplo, ElMatrix_z A)

ElError ElHPSDCholeskyDist_s(ElUpperOrLower uplo, ElDistMatrix_s A)

ElError ElHPSDCholeskyDist_d(ElUpperOrLower uplo, ElDistMatrix_d A)

ElError ElHPSDCholeskyDist_c(ElUpperOrLower uplo, ElDistMatrix_c A)

ElError ElHPSDCholeskyDist_z(ElUpperOrLower uplo, ElDistMatrix_z A)