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.

Main header file

Subroutines

Test driver

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, DistMatrix<int, UPerm, STAR> &p)

Performs Cholesky factorization with full (diagonal) pivoting. On exit, the vector \(p\) consists of the nonzero column indices of the permutation matrix \(P\) such that \(P A P^T = L L^H = U^H U\).

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

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

Updates the 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.

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

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

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

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

Solve linear systems using an unpivoted Cholesky factorization.

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

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

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]

Main header file

Subroutines

Test driver

Example driver

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, UPerm, STAR> &p, LDLPivotType pivotType = BUNCH_KAUFMAN_A)

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

Returns a pivoted LDL factorization, where the vector \(p\) contains the column indices of the nonzero entries of the permutation matrix \(P\) such that \(PAP^T\) equals either \(LDL^T\) or \(LDL^H\), where \(D\) is quasi-diagonal. 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

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

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

Solve linear systems using an unpivoted LDL factorization.

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

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

Solve linear systems using a pivoted LDL factorization.

LU factorization

Main header file

Subroutines

Test driver

Example driver

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, UPerm, STAR> &p)

Overwrites the matrix \(A\) with the LU decomposition of \(PA\), where \(P\) is represented by the permutation vector p, which consists of the column indices of the nonzero entry in each row of \(P\).

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

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

Overwrites the matrix \(A\) with the LU decomposition of \(PAQ^T\), where \(P\) and \(Q\) are represented by the permutation vectors p and q, which consist of the column indices of the nonzero entry in each row of \(P\) and \(Q\), respectively.

void LUMod(Matrix<F> &A, Matrix<int> &p, const Matrix<F> &u, const Matrix<F> &v, bool conjugate = true, Base<F> tau = 0.1)

void LUMod(DistMatrix<F> &A, DistMatrix<int, UPerm, STAR> &p, const DistMatrix<F> &u, const DistMatrix<F> &v, bool conjugate = true, Base<F> tau = 0.1)

Modify an existing LU factorization, \(A = P^T L U\), to incorporate the rank-one update \(A + u v^T\) or \(A + u v^H\). This algorithm only requires a quadratic number of operations.

Note

The current implementation has only been tested for square matrices.

lu namespace

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

void lu::SolveAfter(Orientation orientation, const DistMatrix<F> &A, DistMatrix<F> &B)

Update \(B := A^{-1} B\), \(B := A^{-T} B\), or \(B := A^{-H} B\), where \(A\) has been overwritten with its LU factors (without partial pivoting).

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

void lu::SolveAfter(Orientation orientation, const DistMatrix<F> &A, const DistMatrix<int, UPerm, STAR> &p, DistMatrix<F> &B)

HERE Update \(B := A^{-1} B\), \(B := A^{-T} B\), or \(B := A^{-H} B\), where \(A\) has been overwritten with its LU factors with partial pivoting, which satisfy \(P A = L U\), where the permutation matrix \(P\) is represented by the pivot vector p.

void lu::SolveAfter(Orientation orientation, const Matrix<F> &A, const Matrix<int> &p, const Matrix<int> &q, Matrix<F> &B)

void lu::SolveAfter(Orientation orientation, const DistMatrix<F> &A, const DistMatrix<int, UPerm, STAR> &p, const DistMatrix<int, UPerm, STAR> &q, DistMatrix<F> &B)

Update \(B := A^{-1} B\), \(B := A^{-T} B\), or \(B := A^{-H} B\), where \(A\) has been overwritten with its LU factors with full pivoting, which satisfy \(P A Q = L U\), where the permutation matrices \(P\) and \(Q\) are represented by the pivot vector p and q, respectively.

LQ factorization

Main header file

Subroutines

Test driver

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 (and preceding unitary diagonal matrix forcing \(L\) to have a positive diagonal, defined by the vector d) 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, Matrix<Base<F>> &d)

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

Overwrite the matrix \(A\) with \(L\) and the Householder reflectors representing \(\hat Q\). The scalings for the Householder reflectors are stored in the vector t and the diagonal matrix which forces \(L\) to be positive in d.

lq namespace

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

void lq::ApplyQ(LeftOrRight side, Orientation orientation, const DistMatrix<F> &A, const DistMatrix<F, Ut, Vt> &t, const DistMatrix<Base<F>, Ud, Vd> &d, DistMatrix<F> &B)

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

void lq::SolveAfter(Orientation orientation, const Matrix<F> &A, const Matrix<F> &t, const Matrix<Base<F>> &d, const Matrix<F> &B, Matrix<F> &X)

void lq::SolveAfter(Orientation orientation, const DistMatrix<F> &A, const DistMatrix<F, MD, STAR> &t, const DistMatrix<Base<F>, MD, STAR> &d, const DistMatrix<F> &B, DistMatrix<F> &X)

Solves a set of linear systems using an existing packed LQ factorization given by \(A\) and the vectors \(t\) and \(d\). \(B\) is the matrix of input vectors and \(X\) is the matrix of solutions.

QR factorization

Main header file

Subroutines

Test driver

Example driver

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 (this unitary matrix is scaled from the right by a unitary diagonal matrix with entries given by d so that \(R\) has a positive diagonal).

void QR(Matrix<F> &A)

void QR(DistMatrix<F> &A)

Overwrite \(A\) with \(R\).

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

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

Overwrite the matrix \(A\) with both \(R\) and the Householder reflectors (and subsequent unitary diagonal matrix defined by the vector, d) 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, UPerm, STAR> &p)

Overwrite \(A\) with the \(R\) from a column-pivoted QR factorization, \(A P = Q R\). The permutation matrix \(P\) is represented via the permutation vector \(p\), which contains the column indices of the nonzero entry in each row of \(P\).

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

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

Overwrite \(A\) with both the \(R\) and (scaled) Householder reflectors from a column-pivoted QR factorization.

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::Explicit(Matrix<F> &A, Matrix<F> &R, Matrix<Int> &p)

void qr::Explicit(DistMatrix<F> &A, DistMatrix<F> &R, DistMatrix<int, UPerm, STAR> &p)

Return representations of all matrices of the pivoted QR factorization (note that the pivot vector is returned, not the full pivot matrix).

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

void qr::ApplyQ(LeftOrRight side, Orientation orientation, const DistMatrix<F> &A, const DistMatrix<F, Ut, Vt> &t, const DistMatrix<Base<F>, Ud, Vd> &d, DistMatrix<F> &B)

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

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

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

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

int qr::BusingerGolub(DistMatrix<F> &A, DistMatrix<F, MD, STAR> &t, DistMatrix<Base<F>, MD, STAR> &d, DistMatrix<int, UPerm, 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. The return value is the number of performed iterations.

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

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

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

int qr::BusingerGolub(DistMatrix<F> &A, DistMatrix<F, MD, STAR> &t, DistMatrix<Base<F>, MD, STAR> &d, DistMatrix<int, UPerm, STAR> &p, int numSteps)

Same as above, but only execute a fixed number of steps of the rank-revealing factorization. The return value is the number of performed iterations.

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

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

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

int qr::BusingerGolub(DistMatrix<F> &A, DistMatrix<F, MD, STAR> &t, DistMatrix<int, UPerm, 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. The return value is the number of performed iterations.

void qr::SolveAfter(Orientation orientation, const Matrix<F> &A, const Matrix<F> &t, const Matrix<Base<F>> &d, const Matrix<F> &B, Matrix<F> &X)

void qr::SolveAfter(Orientation orientation, const DistMatrix<F> &A, const DistMatrix<F, MD, STAR> &t, const DistMatrix<Base<F>, MD, STAR> &d, const DistMatrix<F> &B, DistMatrix<F> &X)

Solves a set of linear systems using an existing packed QR factorization given by \(A\) and the vectors \(t\) and \(d\). \(B\) is the matrix of input vectors and \(X\) is the matrix of solutions.

template<>
type TreeData<F>

Matrix<F> QR0

Initial QR factorization

Matrix<F> t0

Phases from initial QR factorization

Matrix<Base<F>> d0

Signature (-1,+1) which scales the Householder matrix from the right.

std::vector<Matrix<F>> QRList

Factorizations within reduction tree

std::vector<Matrix<F>> tList

Phases within reduction tree

std::vector<Matrix<Base<F>>> dList

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

Generalized QR factorization

Implementation

The generalized QR factorization of a pair of matrices \((A,B)\) is analogous to a QR factorization of \(B^{-1} A\) but does not require that \(B\) is square or invertible: unitary matrices \(Q\) and \(Z\), and (right) upper-triangular matrices \(R\) and \(T\), are computed such that

\[A = Q R\]

and

\[B = Q T Z.\]

Thus, if \(B\) was square and invertible, then the QR factorization of \(B^{-1} A\) would be given by \(Z^H (T^{-1} R)\).

void GQR(Matrix<F> &A, Matrix<F> &B)

void GQR(DistMatrix<F> &A, DistMatrix<F> &B)

Overwrite \(A\) with \(R\) and \(B\) with \(T\).

void GQR(Matrix<F> &A, Matrix<F> &tA, Matrix<Base<F>> &dA, Matrix<F> &B, Matrix<F> &tB, Matrix<Base<F>> &dB)

void GQR(DistMatrix<F> &A, DistMatrix<F, MD, STAR> &tA, DistMatrix<Base<F>, MD, STAR> &dA, DistMatrix<F> &B, DistMatrix<F, MD, STAR> &tB, DistMatrix<Base<F>, MD, STAR> &dB)

Overwrite \(A\) with both \(R\) and the (scaled) Householder vectors which, along with the scalings \(tA\) and sign changes \(dA\), define \(Q\). Likewise, \(B\) is overwritten with both \(T\) and the (scaled) Householder vectors which define \(Z\).

RQ factorization

Main header file

Subroutines

Test driver

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, Matrix<Base<F>> &d)

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

Overwrite the matrix \(A\) with \(R\) and the Householder reflectors representing \(\hat Q\). The scalings for the Householder reflectors are stored in the vector t and the unitary diagonal matrix which forces \(R\) to be positive is defined by the vector d.

rq namespace

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

void rq::ApplyQ(LeftOrRight side, Orientation orientation, const DistMatrix<F> &A, const DistMatrix<F, Ut, Vt> &t, const DistMatrix<Base<F>, Ud, Vd> &d, DistMatrix<F> &B)

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

void rq::SolveAfter(Orientation orientation, const Matrix<F> &A, const Matrix<F> &t, const Matrix<Base<F>> &d, const Matrix<F> &B, Matrix<F> &X)

void rq::SolveAfter(Orientation orientation, const DistMatrix<F> &A, const DistMatrix<F, MD, STAR> &t, const DistMatrix<Base<F>, MD, STAR> &d, const DistMatrix<F> &B, DistMatrix<F> &X)

Solves a set of linear systems using an existing packed RQ factorization given by \(A\) and the vectors \(t\) and \(d\). \(B\) is the matrix of input vectors and \(X\) is the matrix of solutions.

Generalized RQ factorization

Implementation

The generalized RQ factorization of a pair of matrices \((A,B)\) is analogous to an RQ factorization of \(A B^{-1}\) but does not require that \(B\) is square or invertible: unitary matrices \(Q\) and \(Z\), and (right) upper-triangular matrices \(R\) and \(T\), are computed such that

\[A = R Q\]

and

\[B = Z T Q.\]

Thus, is \(B\) was square and invertible, then the RQ factorization of \(A B^{-1}\) would be given by \((R T^{-1}) Z^H\).

void GRQ(Matrix<F> &A, Matrix<F> &B)

void GRQ(DistMatrix<F> &A, DistMatrix<F> &B)

Overwrite \(A\) with \(R\) and \(B\) with \(T\).

void GRQ(Matrix<F> &A, Matrix<F> &tA, Matrix<Base<F>> &dA, Matrix<F> &B, Matrix<F> &tB, Matrix<Base<F>> &dB)

void GRQ(DistMatrix<F> &A, DistMatrix<F, MD, STAR> &tA, DistMatrix<Base<F>, MD, STAR> &dA, DistMatrix<F> &B, DistMatrix<F, MD, STAR> &tB, DistMatrix<Base<F>, MD, STAR> &dB)

Overwrite \(A\) with both \(R\) and the (scaled) Householder vectors which, along with the scalings \(tA\) and sign changes \(dA\), define \(Q\). Likewise, \(B\) is overwritten with both \(T\) and the (scaled) Householder vectors which define \(Z\).

Interpolative Decomposition (ID)

Implementation

Example driver

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, UPerm, 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, UPerm, 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

Implementation

Example driver

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, UPerm, STAR> &pR, DistMatrix<int, UPerm, 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.