Matrix properties

Condition number

The condition number of a matrix with respect to a particular norm is

\[\kappa(A) = \|A\| \|A^{-1}\|,\]

with the most common choice being the matrix two-norm.

Base<F> Condition(const Matrix<F> &A, NormType type = TWO_NORM)
Base<F> Condition(const DistMatrix<F, U, V> &A, NormType type = TWO_NORM)

Returns the condition number with respect to the specified norm (one, two, or Frobenius).

Base<F> FrobeniusCondition(const Matrix<F> &A)
Base<F> FrobeniusCondition(const DistMatrix<F, U, V> &A)

Returns the condition number with respect to the Frobenius norm.

Base<F> InfinityCondition(const Matrix<F> &A)
Base<F> InfinityCondition(const DistMatrix<F, U, V> &A)

Returns the condition number with respect to the infinity norm.

Base<F> MaxCondition(const Matrix<F> &A)
Base<F> MaxCondition(const DistMatrix<F, U, V> &A)

Returns the condition number with respect to the entrywise maximum norm.

Base<F> OneCondition(const Matrix<F> &A)
Base<F> OneCondition(const DistMatrix<F, U, V> &A)

Returns the condition number with respect to the one norm.

Base<F> TwoCondition(const Matrix<F> &A)
Base<F> TwoCondition(const DistMatrix<F, U, V> &A)

Returns the condition number with respect to the two norm.

Determinant

Though there are many different possible definitions of the determinant of a matrix \(A \in \mathbb{F}^{n \times n}\), the simplest one is in terms of the product of the eigenvalues (including multiplicity):

\[\mbox{det}(A) = \prod_{i=0}^{n-1} \lambda_i.\]

General

Since \(\mbox{det}(AB)=\mbox{det}(A)\mbox{det}(B)\), we can compute the determinant of an arbitrary matrix in \(\mathcal{O}(n^3)\) work by computing its LU decomposition (with partial pivoting), \(PA=LU\), recognizing that \(\mbox{det}(P)=\pm 1\) (the signature of the permutation), and computing

\[\mbox{det}(A) = \mbox{det}(P)\mbox{det}(L)\mbox{det}(U) = \mbox{det}(P) \prod_{i=0}^{n-1} \upsilon_{i,i} = \pm \prod_{i=0}^{n-1} \upsilon_{i,i},\]

where \(\upsilon_{i,i}\) is the i’th diagonal entry of \(U\).

F Determinant(const Matrix<F> &A)
F Determinant(const DistMatrix<F> &A)
F Determinant(Matrix<F> &A, bool canOverwrite = false)
F Determinant(DistMatrix<F> &A, bool canOverwrite = false)

The determinant of the (fully populated) square matrix A. Some of the variants allow for overwriting the input matrix in order to avoid forming another temporary matrix.

template<>
class SafeProduct<F>

Represents the product of n values as \(\rho \exp(\kappa n)\), where \(|\rho| \le 1\) and \(\kappa \in \mathbb{R}\).

F rho

For nonzero values, rho is the modulus and lies on the unit circle; in order to represent a value that is precisely zero, rho is set to zero.

Base<F> kappa

kappa can be an arbitrary real number.

int n

The number of values in the product.

SafeProduct<F> SafeDeterminant(const Matrix<F> &A)
SafeProduct<F> SafeDeterminant(const DistMatrix<F> &A)
SafeProduct<F> SafeDeterminant(Matrix<F> &A, bool canOverwrite = false)
SafeProduct<F> SafeDeterminant(DistMatrix<F> &A, bool canOverwrite = false)

The determinant of the square matrix A in an expanded form which is less likely to over/under-flow.

HPD

A version of the above determinant specialized for Hermitian positive-definite matrices (which will therefore have all positive eigenvalues and a positive determinant).

Base<F> HPDDeterminant(UpperOrLower uplo, const Matrix<F> &A)
Base<F> HPDDeterminant(UpperOrLower uplo, const DistMatrix<F> &A)
Base<F> HPDDeterminant(UpperOrLower uplo, Matrix<F> &A, bool canOverwrite = false)
Base<F> HPDDeterminant(UpperOrLower uplo, DistMatrix<F> &A, bool canOverwrite = false)

The determinant of the (fully populated) Hermitian positive-definite matrix A. Some of the variants allow for overwriting the input matrix in order to avoid forming another temporary matrix.

SafeProduct<F> SafeHPDDeterminant(UpperOrLower uplo, const Matrix<F> &A)
SafeProduct<F> SafeHPDDeterminant(UpperOrLower uplo, const DistMatrix<F> &A)
SafeProduct<F> SafeHPDDeterminant(UpperOrLower uplo, Matrix<F> &A, bool canOverwrite = false)
SafeProduct<F> SafeHPDDeterminant(UpperOrLower uplo, DistMatrix<F> &A, bool canOverwrite = false)

The determinant of the Hermitian positive-definite matrix A in an expanded form which is less likely to over/under-flow.

Matrix norms

The following routines can return either \(\|A\|_1\), \(\|A\|_\infty\), \(\|A\|_F\) (the Frobenius norm), the maximum entrywise norm, \(\|A\|_2\), or \(\|A\|_*\) (the nuclear/trace norm) of fully-populated matrices.

Base<F> Norm(const Matrix<F> &A, NormType type = FROBENIUS_NORM)
Base<F> Norm(const DistMatrix<F, U, V> &A, NormType type = FROBENIUS_NORM)
Base<F> HermitianNorm(UpperOrLower uplo, const Matrix<F> &A, NormType type = FROBENIUS_NORM)
Base<F> HermitianNorm(UpperOrLower uplo, const DistMatrix<F> &A, NormType type = FROBENIUS_NORM)
Base<F> SymmetricNorm(UpperOrLower uplo, const Matrix<F> &A, NormType type = FROBENIUS_NORM)
Base<F> SymmetricNorm(UpperOrLower uplo, const DistMatrix<F> &A, NormType type = FROBENIUS_NORM)

Compute a norm of a fully-populated or implicitly symmetric/Hermitian (with the data stored in the specified triangle) matrix.

Note

While Norm() supports every type of matrix distribution, HermitianNorm() and SymmetricNorm() currently only support the standard matrix distribution.

Alternatively, one may directly call the following routines (note that the entrywise, KyFan, and Schatten norms have an extra parameter and must be called directly).

Base<F> EntrywiseNorm(const Matrix<F> &A, Base<F> p)
Base<F> EntrywiseNorm(const DistMatrix<F, U, V> &A, Base<F> p)
Base<F> HermitianEntrywiseNorm(UpperOrLower uplo, const Matrix<F> &A, Base<F> p)
Base<F> HermitianEntrywiseNorm(UpperOrLower uplo, const DistMatrix<F> &A, Base<F> p)
Base<F> SymmetricEntrywiseNorm(UpperOrLower uplo, const Matrix<F> &A, Base<F> p)
Base<F> SymmetricEntrywiseNorm(UpperOrLower uplo, const DistMatrix<F> &A, Base<F> p)

The \(\ell_p\) norm of the columns of A stacked into a single vector. Note that the Frobenius norm corresponds to the \(p=2\) case.

Base<F> EntrywiseOneNorm(const Matrix<F> &A)
Base<F> EntrywiseOneNorm(const DistMatrix<F, U, V> &A)
Base<F> HermitianEntrywiseOneNorm(UpperOrLower uplo, const Matrix<F> &A)
Base<F> HermitianEntrywiseOneNorm(UpperOrLower uplo, const DistMatrix<F> &A)
Base<F> SymmetricEntrywiseOneNorm(UpperOrLower uplo, const Matrix<F> &A)
Base<F> SymmetricEntrywiseOneNorm(UpperOrLower uplo, const DistMatrix<F> &A)

The \(\ell_1\) norm of the columns of A stacked into a single vector.

Base<F> FrobeniusNorm(const Matrix<F> &A)
Base<F> FrobeniusNorm(const DistMatrix<F, U, V> &A)
Base<F> HermitianFrobeniusNorm(UpperOrLower uplo, const Matrix<F> &A)
Base<F> HermitianFrobeniusNorm(UpperOrLower uplo, const DistMatrix<F> &A)
Base<F> SymmetricFrobeniusNorm(UpperOrLower uplo, const Matrix<F> &A)
Base<F> SymmetricFrobeniusNorm(UpperOrLower uplo, const DistMatrix<F> &A)

The \(\ell_2\) norm of the singular values (the Schatten norm with \(p=2\)), which can be cheaply computed as the \(\ell_2\) norm of \(\text{vec}(A)\).

Base<F> KyFanNorm(const Matrix<F> &A, int k)
Base<F> KyFanNorm(const DistMatrix<F, U, V> &A, int k)
Base<F> HermitianKyFanNorm(UpperOrLower uplo, const Matrix<F> &A, int k)
Base<F> HermitianKyFanNorm(UpperOrLower uplo, const DistMatrix<F, U, V> &A, int k)
Base<F> SymmetricKyFanNorm(UpperOrLower uplo, const Matrix<F> &A, int k)
Base<F> SymmetricKyFanNorm(UpperOrLower uplo, const DistMatrix<F, U, V> &A, int k)

The sum of the largest k singular values.

Base<F> InfinityNorm(const Matrix<F> &A)
Base<F> InfinityNorm(const DistMatrix<F, U, V> &A)
Base<F> HermitianInfinityNorm(UpperOrLower uplo, const Matrix<F> &A)
Base<F> HermitianInfinityNorm(UpperOrLower uplo, const DistMatrix<F> &A)
Base<F> SymmetricInfinityNorm(UpperOrLower uplo, const Matrix<F> &A)
Base<F> SymmetricInfinityNorm(UpperOrLower uplo, const DistMatrix<F> &A)

The maximum \(\ell_1\) norm of the rows of A. In the symmetric and Hermitian cases, this is equivalent to the \(\|\cdot \|_1\) norm.

Base<F> MaxNorm(const Matrix<F> &A)
Base<F> MaxNorm(const DistMatrix<F, U, V> &A)
Base<F> HermitianMaxNorm(UpperOrLower uplo, const Matrix<F> &A)
Base<F> HermitianMaxNorm(UpperOrLower uplo, const DistMatrix<F> &A)
Base<F> SymmetricMaxNorm(UpperOrLower uplo, const Matrix<F> &A)
Base<F> SymmetricMaxNorm(UpperOrLower uplo, const DistMatrix<F> &A)

The maximum absolute value of the matrix entries.

Base<F> NuclearNorm(const Matrix<F> &A)
Base<F> NuclearNorm(const DistMatrix<F, U, V> &A)
Base<F> HermitianNuclearNorm(UpperOrLower uplo, const Matrix<F> &A)
Base<F> HermitianNuclearNorm(UpperOrLower uplo, const DistMatrix<F, U, V> &A)
Base<F> SymmetricNuclearNorm(UpperOrLower uplo, const Matrix<F> &A)
Base<F> SymmetricNuclearNorm(UpperOrLower uplo, const DistMatrix<F, U, V> &A)

The sum of the singular values. This is equivalent to both the KyFan norm with \(k=n\) and the Schatten norm with \(p=1\). Note that the nuclear norm is dual to the two-norm, which is the Schatten norm with \(p=\infty\).

Base<F> OneNorm(const Matrix<F> &A)
Base<F> OneNorm(const DistMatrix<F, U, V> &A)
Base<F> HermitianOneNorm(UpperOrLower uplo, const Matrix<F> &A)
Base<F> HermitianOneNorm(UpperOrLower uplo, const DistMatrix<F> &A)
Base<F> SymmetricOneNorm(UpperOrLower uplo, const Matrix<F> &A)
Base<F> SymmetricOneNorm(UpperOrLower uplo, const DistMatrix<F> &A)

The maximum \(\ell_1\) norm of the columns of A. In the symmetric and Hermitian cases, this is equivalent to the \(\| \cdot \|_\infty\) norm.

Base<F> SchattenNorm(const Matrix<F> &A, Base<F> p)
Base<F> SchattenNorm(const DistMatrix<F, U, V> &A, Base<F> p)
Base<F> HermitianSchattenNorm(UpperOrLower uplo, const Matrix<F> &A, Base<F> p)
Base<F> HermitianSchattenNorm(UpperOrLower uplo, const DistMatrix<F, U, V> &A, Base<F> p)
Base<F> SymmetricSchattenNorm(UpperOrLower uplo, const Matrix<F> &A, Base<F> p)
Base<F> SymmetricSchattenNorm(UpperOrLower uplo, const DistMatrix<F, U, V> &A, Base<F> p)

The \(\ell_p\) norm of the singular values.

Base<F> TwoNorm(const Matrix<F> &A)
Base<F> TwoNorm(const DistMatrix<F, U, V> &A)
Base<F> HermitianTwoNorm(UpperOrLower uplo, const Matrix<F> &A)
Base<F> HermitianTwoNorm(UpperOrLower uplo, const DistMatrix<F, U, V> &A)
Base<F> SymmetricTwoNorm(UpperOrLower uplo, const Matrix<F> &A)
Base<F> SymmetricTwoNorm(UpperOrLower uplo, const DistMatrix<F, U, V> &A)

The maximum singular value. This is equivalent to the KyFan norm with k equal to one and the Schatten norm with \(p=\infty\).

int ZeroNorm(const Matrix<F> &A)
int ZeroNorm(const DistMatrix<F> &A)
int HermitianZeroNorm(const Matrix<F> &A)
int HermitianZeroNorm(const DistMatrix<F> &A)
int SymmetricZeroNorm(const Matrix<F> &A)
int SymmetricZeroNorm(const DistMatrix<F> &A)

Return the number of nonzero entries in the matrix.

Two-norm estimates

Base<F> TwoNormEstimate(Matrix<F> &A, Base<F> tol = 1e-6)
Base<F> TwoNormEstimate(DistMatrix<F> &A, Base<F> tol = 1e-6)
Base<F> HermitianTwoNormEstimate(Matrix<F> &A, Base<F> tol = 1e-6)
Base<F> HermitianTwoNormEstimate(DistMatrix<F> &A, Base<F> tol = 1e-6)
Base<F> SymmetricTwoNormEstimate(Matrix<F> &A, Base<F> tol = 1e-6)
Base<F> SymmetricTwoNormEstimate(DistMatrix<F> &A, Base<F> tol = 1e-6)

Return an estimate for the two-norm which should be accurate within a factor of \(n\) times the specified tolerance.

Pseudospectra

The \(\epsilon\)-pseudospectrum of a square matrix \(A\) is the set of all shifts \(z\) such that \(\hat A - z\) is singular for some \(\hat A\) such that \(\| \hat A - A \|_2 < \epsilon\). In other words, \(z\) is in the \(\epsilon\)-pseudospectrum of \(A\) if the smallest singular value of \(A - z\) is less than \(\epsilon\).

The method used by Elemental is a high-performance improvement upon the triangularization followed by inverse-iteration approach suggested by Shiu-Hong Lui in Computation of pseudospectra by continuation (please see Trefethen’s Computation of pseudospectra for a comprehensive review). In particular, Elemental begins by computing the Schur decomposition of the given matrix, which preserves the \(\epsilon\)-pseudospectrum, up to round-off error, and then simultaneously performs many Lanczos decompositions on the inverse normal matrix for each shift in a manner which communicates no more data than a standard triangular solve with many right-hand sides. Converged pseudospectrum estimates are deflated after convergence.

Matrix<int> Pseudospectrum(const Matrix<F> &A, const Matrix<Complex<Base<F>>> &shifts, Matrix<Base<F>> &invNorms, bool lanczos = true, bool deflate = true, int maxIts = 1000, Base<F> tol = 1e-6, bool progress = false)
DistMatrix<int, VR, STAR> Pseudospectrum(const DistMatrix<F> &A, const DistMatrix<Complex<Base<F>>, VR, STAR> &shifts, DistMatrix<Base<F>, VR, STAR> &invNorms, bool lanczos = true, bool deflate = true, int maxIts = 1000, Base<F> tol = 1e-6, bool progress = false)
Matrix<int> TriangularPseudospectrum(const Matrix<F> &A, const Matrix<Complex<Base<F>>> &shifts, Matrix<Base<F>> &invNorms, bool lanczos = true, bool deflate = true, int maxIts = 1000, Base<F> tol = 1e-6, bool progress = false)
DistMatrix<int, VR, STAR> TriangularPseudospectrum(const DistMatrix<F> &A, const DistMatrix<Complex<Base<F>>, VR, STAR> &shifts, DistMatrix<Base<F>, VR, STAR> &invNorms, bool lanczos = true, bool deflate = true, int maxIts = 1000, Base<F> tol = 1e-6, bool progress = false)

Returns the norms of the shifted inverses in the vector invNorms for a given set of shifts. The returned integer vector is a list of the number of iterations required for convergence of each shift.

Matrix<int> Pseudospectrum(const Matrix<F> &A, Matrix<Base<F>> &invNormMap, Complex<Base<F>> center, int xSize, int ySize, bool lanczos = true, bool deflate = true, int maxIts = 1000, Base<F> tol = 1e-6, bool progress = false)
DistMatrix<int> Pseudospectrum(const DistMatrix<F> &A, DistMatrix<Base<F>> &invNormMap, Complex<Base<F>> center, int xSize, int ySize, bool lanczos = true, bool deflate = true, int maxIts = 1000, Base<F> tol = 1e-6, bool progress = false)
Matrix<int> TriangularPseudospectrum(const Matrix<F> &A, Matrix<Base<F>> &invNormMap, Complex<Base<F>> center, int xSize, int ySize, bool lanczos = true, bool deflate = true, int maxIts = 1000, Base<F> tol = 1e-6, bool progress = false)
DistMatrix<int> TriangularPseudospectrum(const DistMatrix<F> &A, DistMatrix<Base<F>> &invNormMap, Complex<Base<F>> center, int xSize, int ySize, bool lanczos = true, bool deflate = true, int maxIts = 1000, Base<F> tol = 1e-6, bool progress = false)

Returns the norms of the shifted inverses over a 2D grid (in the matrix invNormMap) with the specified x and y resolutions. The width of the grid in the complex plane is determined based upon the one and two norms of the Schur factor. The returned integer matrix corresponds to the number of iterations required for convergence at each shift in the 2D grid.

Matrix<int> Pseudospectrum(const Matrix<F> &A, Matrix<Base<F>> &invNormMap, Complex<Base<F>> center, Base<F> xWidth, Base<F> yWidth, int xSize, int ySize, bool lanczos = true, bool deflate = true, int maxIts = 1000, Base<F> tol = 1e-6, bool progress = false)
DistMatrix<int> Pseudospectrum(const DistMatrix<F> &A, DistMatrix<Base<F>> &invNormMap, Complex<Base<F>> center, Base<F> xWidth, Base<F> yWidth, int xSize, int ySize, bool lanczos = true, bool deflate = true, int maxIts = 1000, Base<F> tol = 1e-6, bool progress = false)
Matrix<int> TriangularPseudospectrum(const Matrix<F> &A, Matrix<Base<F>> &invNormMap, Complex<Base<F>> center, Base<F> xWidth, Base<F> yWidth, int xSize, int ySize, bool lanczos = true, bool deflate = true, int maxIts = 1000, Base<F> tol = 1e-6, bool progress = false)
DistMatrix<int> TriangularPseudospectrum(const DistMatrix<F> &A, DistMatrix<Base<F>> &invNormMap, Complex<Base<F>> center, Base<F> xWidth, Base<F> yWidth, int xSize, int ySize, bool lanczos = true, bool deflate = true, int maxIts = 1000, Base<F> tol = 1e-6, bool progress = false)

Same as above, but the x and y widths of the 2D grid in the complex plane are manually specified.

Trace

The two equally useful definitions of the trace of a square matrix \(A \in \mathbb{F}^{n \times n}\) are

\[\mbox{tr}(A) = \sum_{i=0}^{n-1} \alpha_{i,i} = \sum_{i=0}^{n-1} \lambda_i,\]

where \(\alpha_{i,i}\) is the i’th diagonal entry of \(A\) and \(\lambda_i\) is the i’th eigenvalue (counting multiplicity) of \(A\).

Clearly the former equation is easier to compute, but the latter is an important characterization.

F Trace(const Matrix<F> &A)
F Trace(const DistMatrix<F> &A)

Return the trace of the square matrix A.

HermitianInertia

void HermitianInertia(UpperOrLower uplo, Matrix<F> &A, LDLPivotType pivotType = BUNCH_PARLETT)
void HermitianInertia(UpperOrLower uplo, DistMatrix<F> &A, LDLPivotType pivotType = BUNCH_PARLETT)

Returns the triplet containing the number of positive, negative, and zero eigenvalues of the Hermitian matrix by analyzing the block diagonal resulting from a pivoted LDL factorization.