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