Matrix functions

Direct inversion

Triangular

Inverts a (possibly unit-diagonal) triangular matrix in-place.

void TriangularInverse(UpperOrLower uplo, UnitOrNonUnit diag, Matrix<F> &A)
void TriangularInverse(UpperOrLower uplo, UnitOrNonUnit diag, DistMatrix<F> &A)

Inverts the triangle of A specified by the parameter uplo; if diag is set to UNIT, then A is treated as unit-diagonal.

General

This routine computes the in-place inverse of a general fully-populated (invertible) matrix \(A\) as

\[A^{-1} = U^{-1} L^{-1} P,\]

where \(PA=LU\) is the result of LU factorization with partial pivoting. The algorithm essentially factors \(A\), inverts \(U\) in place, solves against \(L\) one block column at a time, and then applies the row pivots in reverse order to the columns of the result.

void Inverse(Matrix<F> &A)
void Inverse(DistMatrix<F> &A)

Overwrites the general matrix A with its inverse.

Symmetric/Hermitian

void SymmetricInverse(UpperOrLower uplo, Matrix<F> &A, bool conjugate = false, LDLPivotType pivotType = BUNCH_KAUFMAN_A)
void SymmetricInverse(UpperOrLower uplo, DistMatrix<F> &A, bool conjugate = false, LDLPivotType pivotType = BUNCH_KAUFMAN_A)

Invert a symmetric or Hermitian matrix using a pivoted LDL factorization.

void HermitianInverse(UpperOrLower uplo, Matrix<F> &A, bool conjugate = false, LDLPivotType pivotType = BUNCH_KAUFMAN_A)
void HermitianInverse(UpperOrLower uplo, DistMatrix<F> &A, bool conjugate = false, LDLPivotType pivotType = BUNCH_KAUFMAN_A)

Invert a Hermitian matrix using a pivoted LDL factorization.

HPD

This routine uses a custom algorithm for computing the inverse of a Hermitian positive-definite matrix \(A\) as

\[A^{-1} = (L L^H)^{-1} = L^{-H} L^{-1},\]

where \(L\) is the lower Cholesky factor of \(A\) (the upper Cholesky factor is computed in the case of upper-triangular storage). Rather than performing Cholesky factorization, triangular inversion, and then the Hermitian triangular outer product in sequence, this routine makes use of the single-sweep algorithm described in Bientinesi et al.’s “Families of algorithms related to the inversion of a symmetric positive definite matrix”, in particular, the variant 2 algorithm from Fig. 9.

If the matrix is found to not be HPD, then a NonHPDMatrixException is thrown.

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

Overwrite the uplo triangle of the HPD matrix A with the same triangle of the inverse of A.

Hermitian functions

Reform the matrix with the eigenvalues modified by a user-defined function. When the user-defined function is real-valued, the result will remain Hermitian, but when the function is complex-valued, the result is best characterized as normal.

When the user-defined function, say \(f\), is analytic, we can say much more about the result: if the eigenvalue decomposition of the Hermitian matrix \(A\) is \(A=Z \Lambda Z^H\), then

\[f(A) = f(Z \Lambda Z^H) = Z f(\Lambda) Z^H.\]

Two important special cases are \(f(\lambda) = \exp(\lambda)\) and \(f(\lambda)=\exp(i \lambda)\), where the former results in a Hermitian matrix and the latter in a normal (in fact, unitary) matrix.

Note

Since Elemental currently depends on PMRRR for its tridiagonal eigensolver, only double-precision results are supported as of now.

void RealHermitianFunction(UpperOrLower uplo, Matrix<F> &A, const RealFunctor &f)
void RealHermitianFunction(UpperOrLower uplo, DistMatrix<F> &A, const RealFunctor &f)

Modifies the eigenvalues of the passed-in Hermitian matrix by replacing each eigenvalue \(\lambda_i\) with \(f(\lambda_i) \in \mathbb{R}\). RealFunctor is any class which has the member function Real operator()( Real omega ) const.

void ComplexHermitianFunction(UpperOrLower uplo, Matrix<Complex<Real>> &A, const ComplexFunctor &f)
void ComplexHermitianFunction(UpperOrLower uplo, DistMatrix<Complex<Real>> &A, const ComplexFunctor &f)

Modifies the eigenvalues of the passed-in complex Hermitian matrix by replacing each eigenvalue \(\lambda_i\) with \(f(\lambda_i) \in \mathbb{C}\). ComplexFunctor can be any class which has the member function Complex<Real> operator()( Real omega ) const.

TODO: A version of ComplexHermitianFunction which begins with a real matrix

Pseudoinverse

Pseudoinverse(Matrix<F> &A, Base<F> tolerance = 0)
Pseudoinverse(DistMatrix<F> &A, Base<F> tolerance = 0)

Computes the pseudoinverse of a general matrix through computing its SVD, modifying the singular values with the function

\[\begin{split}f(\sigma) = \left\{\begin{array}{cc} 1/\sigma, & \sigma \ge \epsilon \, n \, \| A \|_2 \\ 0, & \mbox{otherwise} \end{array}\right.,\end{split}\]

where \(\epsilon\) is the relative machine precision, \(n\) is the height of \(A\), and \(\| A \|_2\) is the maximum singular value. If a nonzero value for tolerance was specified, it is used instead of \(\epsilon n \| A \|_2\).

HermitianPseudoinverse(UpperOrLower uplo, Matrix<F> &A, Base<F> tolerance = 0)
HermitianPseudoinverse(UpperOrLower uplo, DistMatrix<F> &A, Base<F> tolerance = 0)

Computes the pseudoinverse of a Hermitian matrix through a customized version of RealHermitianFunction() which used the eigenvalue mapping function

\[\begin{split}f(\omega) = \left\{\begin{array}{cc} 1/\omega, & |\omega| \ge \epsilon \, n \, \| A \|_2 \\ 0, & \mbox{otherwise} \end{array}\right.,\end{split}\]

where \(\epsilon\) is the relative machine precision, \(n\) is the height of \(A\), and \(\| A \|_2\) can be computed as the maximum absolute value of the eigenvalues of \(A\). If a nonzero value for tolerance is specified, it is used instead of \(\epsilon n \| A \|_2\).

Square root

A matrix \(B\) satisfying

\[B^2 = A\]

is referred to as the square-root of the matrix \(A\). Such a matrix is guaranteed to exist as long as \(A\) is diagonalizable: if \(A = X \Lambda X^{-1}\), then we may put

\[B = X \sqrt{\Lambda} X^{-1},\]

where each eigenvalue \(\lambda = r e^{i\theta}\) maps to \(\sqrt{\lambda} = \sqrt{r} e^{i\theta/2}\).

void SquareRoot(Matrix<F> &A)
void SquareRoot(DistMatrix<F> &A)

Currently uses a Newton iteration to compute the general matrix square-root. See square_root::Newton for the more detailed interface.

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

Computes the Hermitian EVD, square-roots the eigenvalues, and then reforms the matrix. If any of the eigenvalues were sufficiently negative, a NonHPSDMatrixException is thrown.

TODO: HermitianSquareRoot

square_root namespace

int square_root::Newton(Matrix<F> &A, int maxIts = 100, Base<F> tol = 0)
int square_root::Newton(DistMatrix<F> &A, int maxIts = 100, Base<F> tol = 0)

Performs at most maxIts Newton steps in an attempt to compute the matrix square-root within the specified tolerance, which defaults to \(n \epsilon\), where \(n\) is the matrix height and \(\epsilon\) is the machine precision.

Sign

The matrix sign function can be written as

\[\text{sgn}(A) = A(A^2)^{-1/2},\]

as long as \(A\) does not have any pure-imaginary eigenvalues.

void Sign(Matrix<F> &A)
void Sign(DistMatrix<F> &A)
void Sign(Matrix<F> &A, Matrix<F> &N)
void Sign(DistMatrix<F> &A, DistMatrix<F> &N)

Compute the matrix sign through a globally-convergent Newton iteration scaled with the Frobenius norm of the iterate and its inverse. Optionally return the full decomposition, \(A=S N\), where \(A\) is overwritten by \(S\).

void HermitianSign(UpperOrLower uplo, Matrix<F> &A)
void HermitianSign(UpperOrLower uplo, DistMatrix<F> &A)
void HermitianSign(UpperOrLower uplo, Matrix<F> &A, Matrix<F> &N)
void HermitianSign(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<F> &N)

Compute the Hermitian EVD, replace the eigenvalues with their sign, and then reform the matrix. Optionally return the full decomposition, \(A=SN\), where \(A\) is overwritten by \(S\). Note that this will also be a polar decomposition.

sign namespace

type sign::Scaling

An enum for specifying the scaling strategy to be used for the Newton iteration for the matrix sign function. It must be either NONE, DETERMINANT, or FROB_NORM (the default).

int sign::Newton(Matrix<F> &A, sign::Scaling scaling = FROB_NORM, int maxIts = 100, Base<F> tol = 0)
int sign::Newton(DistMatrix<F> &A, sign::Scaling scaling = FROB_NORM, int maxIts = 100, Base<F> tol = 0)

Runs a (scaled) Newton iteration for at most maxIts iterations with the specified tolerance, which, if undefined, is set to \(n \epsilon\), where \(n\) is the matrix dimension and \(\epsilon\) is the machine epsilon. The return value is the number of performed iterations.