Factorizations

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.

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

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

Overwrite the uplo triangle of the HPD matrix A with its Cholesky factor.

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

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

Performs Cholesky factorization with full (diagonal) pivoting.

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.

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

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

Overwrite the uplo triangle of the potentially singular matrix A with its Cholesky factor.

cholesky namespace

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

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

Solve linear systems using an unpivoted Cholesky factorization.

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

cholesky::SolveAfter(UpperOrLower uplo, Orientation orientation, const DistMatrix<F> &A, DistMatrix<F> &B, DistMatrix<int, VC, STAR> &p)

Solve linear systems using a pivoted Cholesky factorization.

\(LDL\) factorization

enum LDLPivotType

For specifying a symmetric pivoting strategy. The current (not yet all supported) options include:

enumerator BUNCH_KAUFMAN_A

enumerator BUNCH_KAUFMAN_C

Not yet supported.

enumerator BUNCH_KAUFMAN_D

enumerator BUNCH_KAUFMAN_BOUNDED

Not yet supported.

enumerator BUNCH_PARLETT

class LDLPivot

int nb

int from[2]

void LDLH(Matrix<F> &A, Matrix<F> &dSub, Matrix<int> &p, LDLPivotType pivotType = BUNCH_KAUFMAN_A)

void LDLT(Matrix<F> &A, Matrix<F> &dSub, Matrix<int> &p, LDLPivotType pivotType = BUNCH_KAUFMAN_A)

void LDLH(DistMatrix<F> &A, DistMatrix<F, MD, STAR> &dSub, DistMatrix<int, VC, STAR> &p, LDLPivotType pivotType = BUNCH_KAUFMAN_A)

void LDLT(DistMatrix<F> &A, DistMatrix<F, MD, STAR> &dSub, DistMatrix<int, VC, STAR> &p, LDLPivotType pivotType = BUNCH_KAUFMAN_A)

Pivoted LDL factorization. The Bunch-Kaufman pivoting rules are used within a higher-performance blocked algorithm, whereas the Bunch-Parlett strategy uses an unblocked algorithm.

Though the Cholesky factorization is ideal for most HPD matrices, the unpivoted LDL factorizations exist as slight relaxation of the Cholesky factorization and compute lower-triangular (with unit diagonal) \(L\) and diagonal \(D\) such that \(A = L D L^H\) or \(A = L D L^T\). If a zero pivot is attempted, then a ZeroPivotException will be thrown.

Warning

The following routines do not pivot, so please use with caution.

void LDLH(Matrix<F> &A)

void LDLT(Matrix<F> &A)

void LDLH(DistMatrix<F> &A)

void LDLT(DistMatrix<F> &A)

Overwrite the strictly lower triangle of \(A\) with the strictly lower portion of \(L\) (\(L\) implicitly has ones on its diagonal) and the diagonal with \(D\).

ldl namespace

ldl::SolveAfter(const Matrix<F> &A, Matrix<F> &B, bool conjugated = false)

ldl::SolveAfter(const DistMatrix<F> &A, DistMatrix<F> &B, bool conjugated = false)

Solve linear systems using an unpivoted LDL factorization.

ldl::SolveAfter(const Matrix<F> &A, const Matrix<F> &dSub, const Matrix<int> &p, Matrix<F> &B, bool conjugated = false)

ldl::SolveAfter(const DistMatrix<F> &A, const DistMatrix<F, MD, STAR> &dSub, const DistMatrix<int, VC, STAR> &p, DistMatrix<F> &B, bool conjugated = false)

Solve linear systems using a pivoted LDL factorization.

\(LU\) factorization

Given \(A \in \mathbb{F}^{m \times n}\), an LU factorization (without pivoting) finds a unit lower-trapezoidal \(L \in \mathbb{F}^{m \times \mbox{min}(m,n)}\) and upper-trapezoidal \(U \in \mathbb{F}^{\mbox{min}(m,n) \times n}\) such that \(A=LU\). Since \(L\) is required to have its diaganal entries set to one: the upper portion of \(A\) can be overwritten with U, and the strictly lower portion of \(A\) can be overwritten with the strictly lower portion of \(L\). If \(A\) is found to be numerically singular, then a SingularMatrixException will be thrown.

void LU(Matrix<F> &A)

void LU(DistMatrix<F> &A)

Overwrites \(A\) with its LU decomposition.

Since LU factorization without pivoting is known to be unstable for general matrices, it is standard practice to pivot the rows of \(A\) during the factorization (this is called partial pivoting since the columns are not also pivoted). An LU factorization with partial pivoting therefore computes \(P\), \(L\), and \(U\) such that \(PA=LU\), where \(L\) and \(U\) are as described above and \(P\) is a permutation matrix.

void LU(Matrix<F> &A, Matrix<int> &p)

void LU(DistMatrix<F> &A, DistMatrix<F, VC, STAR> &p)

Overwrites the matrix \(A\) with the LU decomposition of \(PA\), where \(P\) is represented by the pivot vector p.

void LU(Matrix<F> &A, Matrix<int> &p, Matrix<int> &q)

void LU(DistMatrix<F> &A, DistMatrix<F, VC, STAR> &p, DistMatrix<F, VC, STAR> &q)

Overwrites the matrix \(A\) with the LU decomposition of \(PAQ\), where \(P\) is represented by the pivot vector p, and likewise for \(Q\).

\(LQ\) factorization

Given \(A \in \mathbb{F}^{m \times n}\), an LQ factorization typically computes an implicit unitary matrix \(\hat Q \in \mathbb{F}^{n \times n}\) such that \(\hat L \equiv A\hat Q^H\) is lower trapezoidal. One can then form the thin factors \(L \in \mathbb{F}^{m \times \mbox{min}(m,n)}\) and \(Q \in \mathbb{F}^{\mbox{min}(m,n) \times n}\) by setting \(L\) and \(Q\) to first \(\mbox{min}(m,n)\) columns and rows of \(\hat L\) and \(\hat Q\), respectively. Upon completion \(L\) is stored in the lower trapezoid of \(A\) and the Householder reflectors representing \(\hat Q\) are stored within the rows of the strictly upper trapezoid.

void LQ(Matrix<F> &A)

void LQ(DistMatrix<F> &A)

void LQ(Matrix<F> &A, Matrix<F> &t)

void LQ(DistMatrix<F> &A, DistMatrix<F, MD, STAR> &t)

Overwrite the complex matrix \(A\) with \(L\) and the Householder reflectors representing \(\hat Q\). The scalings for the Householder reflectors are stored in the vector t.

lq namespace

void lq::ApplyQ(LeftOrRight side, Orientation orientation, const Matrix<F> &A, const Matrix<F> &t, Matrix<F> &B)

void lq::ApplyQ(LeftOrRight side, Orientation orientation, const DistMatrix<F> &A, const DistMatrix<F, MD, STAR> &t, DistMatrix<F> &B)

void lq::ApplyQ(LeftOrRight side, Orientation orientation, const DistMatrix<F> &A, const DistMatrix<F, STAR, STAR> &t, DistMatrix<F> &B)

Applies the implicitly-defined \(Q\) (or its adjoint) stored within A and t from either the left or the right to \(B\).

\(QR\) factorization

Given \(A \in \mathbb{F}^{m \times n}\), a QR factorization typically computes an implicit unitary matrix \(\hat Q \in \mathbb{F}^{m \times m}\) such that \(\hat R \equiv \hat Q^H A\) is upper trapezoidal. One can then form the thin factors \(Q \in \mathbb{F}^{m \times \mbox{min}(m,n)}\) and \(R \in \mathbb{F}^{\mbox{min}(m,n) \times n}\) by setting \(Q\) and \(R\) to first \(\mbox{min}(m,n)\) columns and rows of \(\hat Q\) and \(\hat R\), respectively. Upon completion \(R\) is stored in the upper trapezoid of \(A\) and the Householder reflectors representing \(\hat Q\) are stored within the columns of the strictly lower trapezoid.

void QR(Matrix<F> &A)

void QR(DistMatrix<F> &A)

void QR(Matrix<F> &A, Matrix<F> &t)

void QR(DistMatrix<F> &A, DistMatrix<F, MD, STAR> &t)

Overwrite the complex matrix \(A\) with \(R\) and the Householder reflectors representing \(\hat Q\). The scalings for the Householder reflectors are stored in the vector t.

void QR(Matrix<F> &A, Matrix<int> &p)

void QR(DistMatrix<F> &A, DistMatrix<int, VR, STAR> &p)

void QR(Matrix<F> &A, Matrix<F> &t, Matrix<int> &p)

void QR(DistMatrix<F> &A, DistMatrix<F, MD, STAR> &t, DistMatrix<int, VR, STAR> &p)

Column-pivoted QR factorization. The current implementation uses Businger-Golub pivoting.

qr namespace

void qr::Explicit(Matrix<F> &A, bool colPiv = false)

void qr::Explicit(DistMatrix<F> &A, bool colPiv = false)

Overwrite \(A\) with the orthogonal matrix from its QR factorization (with or without column pivoting).

void qr::Explicit(Matrix<F> &A, Matrix<F> &R, bool colPiv = false)

void qr::Explicit(DistMatrix<F> &A, DistMatrix<F> &R, bool colPiv = false)

Additionally explicitly return the \(R\) from the QR factorization.

void qr::ApplyQ(LeftOrRight side, Orientation orientation, const Matrix<F> &A, const Matrix<F> &t, Matrix<F> &B)

void qr::ApplyQ(LeftOrRight side, Orientation orientation, const DistMatrix<F> &A, const DistMatrix<F, MD, STAR> &t, DistMatrix<F> &B)

void qr::ApplyQ(LeftOrRight side, Orientation orientation, const DistMatrix<F> &A, const DistMatrix<F, STAR, STAR> &t, DistMatrix<F> &B)

Applies the implicitly-defined \(Q\) (or its adjoint) stored within A and t from either the left or the right to \(B\).

void qr::BusingerGolub(Matrix<F> &A, Matrix<int> &p)

void qr::BusingerGolub(DistMatrix<F> &A, DistMatrix<int, VR, STAR> &p)

void qr::BusingerGolub(Matrix<F> &A, Matrix<F> &t, Matrix<int> &p)

void qr::BusingerGolub(DistMatrix<F> &A, DistMatrix<F, MD, STAR> &t, DistMatrix<int, VR, STAR> &p)

Column-pivoted versions of the above routines which use the Businger/Golub strategy, i.e., the pivot is chosen as the remaining column with maximum two norm.

void qr::BusingerGolub(Matrix<F> &A, Matrix<int> &p, int numSteps)

void qr::BusingerGolub(DistMatrix<F> &A, DistMatrix<int, VR, STAR> &p, int numSteps)

void qr::BusingerGolub(Matrix<F> &A, Matrix<F> &t, Matrix<int> &p, int numSteps)

void qr::BusingerGolub(DistMatrix<F> &A, DistMatrix<F, MD, STAR> &t, DistMatrix<int, VR, STAR> &p, int numSteps)

Same as above, but only execute a fixed number of steps of the rank-revealing factorization.

void qr::BusingerGolub(Matrix<F> &A, Matrix<int> &p, int maxSteps, Base<F> tol)

void qr::BusingerGolub(DistMatrix<F> &A, DistMatrix<int, VR, STAR> &p, int maxSteps, Base<F> tol)

void qr::BusingerGolub(Matrix<F> &A, Matrix<F> &t, Matrix<int> &p, int maxSteps, Base<F> tol)

void qr::BusingerGolub(DistMatrix<F> &A, DistMatrix<F, MD, STAR> &t, DistMatrix<int, VR, STAR> &p, int maxSteps, Base<F> tol)

Either execute maxSteps iterations or stop after the maximum remaining column norm is less than or equal to tol times the maximum original column norm.

template<>
class TreeData<F>

Matrix<F> QR0

Initial QR factorization

Matrix<F> t0

Phases from initial QR factorization

std::vector<Matrix<F>> QRList

Factorizations within reduction tree

std::vector<Matrix<F>> tList

Phases within reduction tree

qr::TreeData<F> qr::TS(const DistMatrix<F, U, STAR> &A)

Forms an implicit tall-skinny QR decomposition.

void qr::ExplicitTS(DistMatrix<F, U, STAR> &A, DistMatrix<F, STAR, STAR> &R)

Forms an explicit QR decomposition using a tall-skinny algorithm: A is overwritten with Q.

qr::ts namespace

DistMatrix<F, STAR, STAR> qr::ts::FormR(const DistMatrix<F, U, STAR> &A, const qr::TreeData<F> &treeData)

Return the R from the QR decomposition.

void qr::ts::FormQ(DistMatrix<F, U, STAR> &A, qr::TreeData<F> &treeData)

Overwrite A with the Q from the QR decomposition.

\(RQ\) factorization

Just like an LQ factorization, but the orthogonalization process starts from the bottom row and produces a much sparser triangular factor when the matrix is wider than it is tall.

void RQ(Matrix<F> &A)

void RQ(DistMatrix<F> &A)

void RQ(Matrix<F> &A, Matrix<F> &t)

void RQ(DistMatrix<F> &A, DistMatrix<F, MD, STAR> &t)

Overwrite the complex matrix \(A\) with \(R\) and the Householder reflectors representing \(\hat Q\). The scalings for the Householder reflectors are stored in the vector t.

rq namespace

void rq::ApplyQ(LeftOrRight side, Orientation orientation, const Matrix<F> &A, const Matrix<F> &t, Matrix<F> &B)

void rq::ApplyQ(LeftOrRight side, Orientation orientation, const DistMatrix<F> &A, const DistMatrix<F, MD, STAR> &t, DistMatrix<F> &B)

void rq::ApplyQ(LeftOrRight side, Orientation orientation, const DistMatrix<F> &A, const DistMatrix<F, STAR, STAR> &t, DistMatrix<F> &B)

Applies the implicitly-defined \(Q\) (or its adjoint) stored within A and t from either the left or the right to \(B\).

Interpolative Decomposition (ID)

Interpolative Decompositions (ID’s) are closely related to pivoted QR factorizations and are useful for representing (approximately) low-rank matrices in terms of linear combinations of a few of their columns, i.e.,

\[A P = \hat{A} \begin{pmatrix} I & Z \end{pmatrix},\]

where \(P\) is a permutation matrix, \(\hat{A}\) is a small set of columns of \(A\), and \(Z\) is an interpolation matrix responsible for representing the remaining columns in terms of the selected columns of \(A\).

void ID(const Matrix<F> &A, Matrix<int> &p, Matrix<F> &Z, int numSteps)

void ID(const DistMatrix<F> &A, DistMatrix<int, VR, STAR> &p, DistMatrix<F, STAR, VR> &Z, int numSteps)

numSteps steps of a pivoted QR factorization are used to return an Interpolative Decomposition of \(A\).

void ID(const Matrix<F> &A, Matrix<int> &p, Matrix<F> &Z, int maxSteps, Base<F> tol)

void ID(const DistMatrix<F> &A, DistMatrix<int, VR, STAR> &p, DistMatrix<F, STAR, VR> &Z, int maxSteps, Base<F> tol)

Either maxSteps steps of a pivoted QR factorization are used, or executation stopped after the maximum remaining column norm was less than or equal to tol times the maximum original column norm.

Skeleton decomposition

Skeleton decompositions are essentially two-sided interpolative decompositions, but the terminology is unfortunately extremely contested. We follow the convention that a skeleton decomposition is an approximation

\[A \approx A_C Z A_R,\]

where \(A_C\) is a (small) selection of columns of \(A\), \(A_R\) is a (small) selection of rows of \(A\), and \(Z\) is a (small) square matrix. When \(Z\) is allowed to be rectangular, it is more common to call this a CUR decomposition.

void Skeleton(const Matrix<F> &A, Matrix<int> &pR, Matrix<int> &pC, Matrix<F> &Z, int maxSteps, Base<F> tol)

void Skeleton(const DistMatrix<F> &A, DistMatrix<int, VR, STAR> &pR, DistMatrix<int, VR, STAR> &pC, int maxSteps, Base<F> tol)

Rather than returning \(A_R\) and \(A_C\), the permutation matrices which implicitly define them are returned instead. At most maxSteps steps of a pivoted QR decomposition will be used in order to generate the row/column subsets, and less steps will be taken if a pivot norm is less than or equal to tolerance times the first pivot norm.