Cholesky factorization

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.

Elemental provides support for versions which perform no pivoting, perform full (in this case, equivalent to diagonal) pivoting, and a version which makes use of a matrix square-root of a Hermitian matrix: 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^H B = (Q R)^H (Q R) = R^H R.\]

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

Python API

Cholesky(uplo, A[, piv=False])
HPSDCholesky(uplo, A)

C++ API

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

No pivoting

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

Full (diagonal) pivoting

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

Use the QR (or LQ) factorization of the matrix square-root of \(A\) to find the Cholesky factor.

C API

Single-precision

ElError ElCholesky_s(ElUpperOrLower uplo, ElMatrix_s A)
ElError ElCholeskyDist_s(ElUpperOrLower uplo, ElDistMatrix_s A)

No pivoting

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

Full (diagonal) pivoting

ElError ElHPSDCholesky_s(ElUpperOrLower uplo, ElMatrix_s A)
ElError ElHPSDCholeskyDist_s(ElUpperOrLower uplo, ElDistMatrix_s A)

Use the QR (or LQ) factorization of the matrix square-root of \(A\) to find the Cholesky factor.

Double-precision

ElError ElCholesky_d(ElUpperOrLower uplo, ElMatrix_d A)
ElError ElCholeskyDist_d(ElUpperOrLower uplo, ElDistMatrix_d A)

No pivoting

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

Full (diagonal) pivoting

ElError ElHPSDCholesky_d(ElUpperOrLower uplo, ElMatrix_d A)
ElError ElHPSDCholeskyDist_d(ElUpperOrLower uplo, ElDistMatrix_d A)

Use the QR (or LQ) factorization of the matrix square-root of \(A\) to find the Cholesky factor.

Single-precision complex

ElError ElCholesky_c(ElUpperOrLower uplo, ElMatrix_c A)
ElError ElCholeskyDist_c(ElUpperOrLower uplo, ElDistMatrix_c A)

No pivoting

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

Full (diagonal) pivoting

ElError ElHPSDCholesky_c(ElUpperOrLower uplo, ElMatrix_c A)
ElError ElHPSDCholeskyDist_c(ElUpperOrLower uplo, ElDistMatrix_c A)

Use the QR (or LQ) factorization of the matrix square-root of \(A\) to find the Cholesky factor.

Double-precision complex

ElError ElCholesky_z(ElUpperOrLower uplo, ElMatrix_z A)
ElError ElCholeskyDist_z(ElUpperOrLower uplo, ElDistMatrix_z A)

No pivoting

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

Full (diagonal) pivoting

ElError ElHPSDCholesky_z(ElUpperOrLower uplo, ElMatrix_z A)
ElError ElHPSDCholeskyDist_z(ElUpperOrLower uplo, ElDistMatrix_z A)

Use the QR (or LQ) factorization of the matrix square-root of \(A\) to find the Cholesky factor.

Solving linear systems with the factorization

After a (possibly pivoted) Cholesky factorization has been formed, it is possible to solve linear systems in quadratic time. The following routines apply the inverse in such a fast manner.

Python API

SolveAfterCholesky(uplo, orient, A, B)

No pivoting

SolveAfterCholesky(uplo, orient, A, p, B)

Full (diagonal) pivoting

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)

No pivoting

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)

Full (diagonal) pivoting

C API

Single-precision

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

No pivoting

ElError ElSolveAfterCholeskyPiv_s(ElUpperOrLower uplo, ElOrientation orientation, ElConstMatrix_s A, ElConstMatrix_i p, ElMatrix_s B)
ElError ElSolveAfterCholeskyPivDist_s(ElUpperOrLower uplo, ElOrientation orientation, ElConstDistMatrix_s A, ElConstDistMatrix_i p, ElDistMatrix_s B)

Full (diagonal) pivoting

Double-precision

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

No pivoting

ElError ElSolveAfterCholeskyPiv_d(ElUpperOrLower uplo, ElOrientation orientation, ElConstMatrix_d A, ElConstMatrix_i p, ElMatrix_d B)
ElError ElSolveAfterCholeskyPivDist_d(ElUpperOrLower uplo, ElOrientation orientation, ElConstDistMatrix_d A, ElConstDistMatrix_i p, ElDistMatrix_d B)

Full (diagonal) pivoting

Single-precision complex

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

No pivoting

ElError ElSolveAfterCholeskyPiv_c(ElUpperOrLower uplo, ElOrientation orientation, ElConstMatrix_c A, ElConstMatrix_i p, ElMatrix_c B)
ElError ElSolveAfterCholeskyPivDist_c(ElUpperOrLower uplo, ElOrientation orientation, ElConstDistMatrix_c A, ElConstDistMatrix_i p, ElDistMatrix_c B)

Full (diagonal) pivoting

Double-precision complex

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

No pivoting

ElError ElSolveAfterCholeskyPiv_z(ElUpperOrLower uplo, ElOrientation orientation, ElConstMatrix_z A, ElConstMatrix_i p, ElMatrix_z B)
ElError ElSolveAfterCholeskyPivDist_z(ElUpperOrLower uplo, ElOrientation orientation, ElConstDistMatrix_z A, ElConstDistMatrix_i p, ElDistMatrix_z B)

Full (diagonal) pivoting

Low-rank updates to a factorization

It is well-known that it is possible to update an existing Cholesky factorization to incorporate a low-rank modification \(\alpha V V^H\) in quadratic time. The following algorithms use Householder transformations for updates (\(\alpha \ge 0\)) and hyperbolic Householder transformations for downdates.

Python API

CholeskyMod(uplo, T, alpha, V)

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

Single-precision

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

Double-precision

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

Single-precision complex

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

Double-precision complex

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