Matrix decompositions

Hermitian eigensolver

Elemental provides a collection of routines for both full and partial solution of the Hermitian eigenvalue problem

\[A Z = Z \Lambda,\]

where A is the given Hermitian matrix, and unitary Z and real diagonal \(\Lambda\) are sought. In particular, with the eigenvalues and corresponding eigenpairs labeled in non-decreasing order, the three basic modes are:

  1. Compute all eigenvalues or eigenpairs, \(\{\lambda_i\}_{i=0}^{n-1}\) or \(\{(z_i,\lambda_i)\}_{i=0}^{n-1}\).
  2. Compute the eigenvalues or eigenpairs with a given range of indices, say \(\{\lambda_i\}_{i=a}^b\) or \(\{(z_i,\lambda_i)\}_{i=a}^b\), with \(0 \le a \le b < n\).
  3. Compute all eigenpairs (or just eigenvalues) with eigenvalues lying in a particular half-open interval, either \(\{\lambda_i \;|\; \lambda_i \in (a,b] \}\) or \(\{ (z_i,\lambda_i) \;|\; \lambda_i \in (a,b] \}\).

As of now, all three approaches start with Householder tridiagonalization (ala HermitianTridiag()) and then call Matthias Petschow and Paolo Bientinesi’s PMRRR for the tridiagonal eigenvalue problem.

Note

Unfortunately, PMRRR currently only supports double-precision problems, and so the parallel versions of these routines are limited to real and complex double-precision matrices.

Note

Please see the Tuning parameters section for information on optimizing the reduction to tridiagonal form, as it is the dominant cost in all of Elemental’s Hermitian eigensolvers.

Full spectrum computation

void HermitianEig(UpperOrLower uplo, Matrix<F> &A, Matrix<Base<F>> &w, SortType sort = UNSORTED)
void HermitianEig(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &w, SortType sort = UNSORTED)

Compute the full set of eigenvalues of the Hermitian matrix A.

void HermitianEig(UpperOrLower uplo, Matrix<F> &A, Matrix<Base<F>> &w, Matrix<F> &Z, SortType sort = UNSORTED)
void HermitianEig(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &w, DistMatrix<F> &Z, SortType sort = UNSORTED)

Compute the full set of eigenpairs of the Hermitian matrix A.

Index-based subset computation

void HermitianEig(UpperOrLower uplo, Matrix<F> &A, Matrix<Base<F>> &w, int a, int b, SortType sort = UNSORTED)
void HermitianEig(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &w, int a, int b, SortType sort = UNSORTED)

Compute the eigenvalues of a Hermitian matrix A with indices in the range \(a,a+1,...,b\).

void HermitianEig(UpperOrLower uplo, Matrix<F> &A, Matrix<Base<F>> &w, Matrix<F> &Z, int a, int b, SortType sort = UNSORTED)
void HermitianEig(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &w, DistMatrix<F> &Z, int a, int b, SortType sort = UNSORTED)

Compute the eigenpairs of a Hermitian matrix A with indices in the range \(a,a+1,...,b\).

Range-based subset computation

void HermitianEig(UpperOrLower uplo, Matrix<F> &A, Matrix<Base<F>> &w, Base<F> a, Base<F> b, SortType sort = UNSORTED)
void HermitianEig(UpperOrLower uplo, DistMatrix<F, STAR, STAR> &A, DistMatrix<Base<F>, STAR, STAR> &w, Base<F> a, Base<F> b, SortType sort = UNSORTED)
void HermitianEig(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &w, Base<F> a, Base<F> b, SortType sort = UNSORTED)

Compute the eigenvalues of a Hermitian matrix A lying in the half-open interval \((a,b]\).

void HermitianEig(UpperOrLower uplo, Matrix<F> &A, Matrix<Base<F>> &w, Matrix<F> &Z, Base<F> a, Base<F> b, SortType sort = UNSORTED)
void HermitianEig(UpperOrLower uplo, DistMatrix<F, STAR, STAR> &A, DistMatrix<Base<F>, STAR, STAR> &w, DistMatrix<F, STAR, STAR> &Z, Base<F> a, Base<F> b, SortType sort = UNSORTED)
void HermitianEig(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &w, DistMatrix<F> &Z, Base<F> a, Base<F> b, SortType sort = UNSORTED)

Compute the eigenpairs of a Hermitian matrix A with eigenvalues lying in the half-open interval \((a,b]\).

Spectral divide and conquer

The primary references for this approach is Demmel et al.’s Fast linear algebra is stable and Nakatsukasa et al.’s Stable and efficient spectral divide and conquer algorithms for the symmetric eigenvalue problem.

void hermitian_eig::SDC(Matrix<F> &A, Matrix<Base<F>> &w, int cutoff = 256, int maxInnerIts = 1, int maxOuterIts = 10, Base<F> relTol = 0)
void hermitian_eig::SDC(DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &w, int cutoff = 256, int maxInnerIts = 1, int maxOuterIts = 10, Base<F> relTol = 0)

Compute the eigenvalues of the matrix \(A\) via a QDWH-based spectral divide and conquer process.

The cutoff controls when the problem is sufficiently small to switch to a standard algorithm, the number of inner iterations is how many attempts to make with the same randomized URV decomposition, and the number of outer iterations is how many random Mobius transformations to try for each spectral split before giving up.

void hermitian_eig::SDC(Matrix<F> &A, Matrix<Base<F>> &w, Matrix<F> &Q, int cutoff = 256, int maxInnerIts = 1, int maxOuterIts = 10, Base<F> relTol = 0)
void hermitian_eig::SDC(DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &w, DistMatrix<F> &Q, int cutoff = 256, int maxInnerIts = 1, int maxOuterIts = 10, Base<F> relTol = 0)

Attempt to also compute the eigenvectors.

Skew-Hermitian eigensolver

Essentially identical to the Hermitian eigensolver, HermitianEig(); for any skew-Hermitian matrix \(G\), \(iG\) is Hermitian, as

\[(iG)^H = -iG^H = iG.\]

This fact implies a fast method for solving skew-Hermitian eigenvalue problems:

  1. Form \(iG\) in \(O(n^2)\) work (switching to complex arithmetic in the real case)
  2. Run a Hermitian eigensolve on \(iG\), yielding \(iG=Z \Lambda Z^H\).
  3. Recognize that \(G=Z (-i \Lambda) Z^H\) provides an EVD of \(G\).

Please see the HermitianEig() documentation for more details.

Note

Unfortunately, PMRRR currently only supports double-precision problems, and so the parallel versions of these routines are limited to real and complex double-precision matrices.

Note

Please see the Tuning parameters section for information on optimizing the reduction to tridiagonal form, as it is the dominant cost in all of Elemental’s Hermitian eigensolvers.

Full spectrum computation

void SkewHermitianEig(UpperOrLower uplo, Matrix<F> &G, Matrix<Base<F>> &wImag, SortType sort = UNSORTED)
void SkewHermitianEig(UpperOrLower uplo, DistMatrix<F> &G, DistMatrix<Base<F>, VR, STAR> &wImag, SortType sort = UNSORTED)

Compute the full set of eigenvalues of the skew-Hermitian matrix G.

void SkewHermitianEig(UpperOrLower uplo, Matrix<F> &G, Matrix<Base<F>> &wImag, Matrix<Complex<Base<F>>> &Z, SortType sort = UNSORTED)
void SkewHermitianEig(UpperOrLower uplo, DistMatrix<F> &G, DistMatrix<Base<F>, VR, STAR> &wImag, DistMatrix<Complex<Base<F>>> &Z, SortType sort = UNSORTED)

Compute the full set of eigenpairs of the skew-Hermitian matrix G.

Index-based subset computation

void SkewHermitianEig(UpperOrLower uplo, Matrix<F> &G, Matrix<Base<F>> &wImag, int a, int b, SortType sort = UNSORTED)
void SkewHermitianEig(UpperOrLower uplo, DistMatrix<F> &G, DistMatrix<Base<F>, VR, STAR> &wImag, int a, int b, SortType sort = UNSORTED)

Compute the eigenvalues of a skew-Hermitian matrix G with indices in the range \(a,a+1,...,b\).

void SkewHermitianEig(UpperOrLower uplo, Matrix<F> &G, Matrix<Base<F>> &wImag, Matrix<Complex<Base<F>>> &Z, int a, int b, SortType sort = UNSORTED)
void SkewHermitianEig(UpperOrLower uplo, DistMatrix<F> &G, DistMatrix<Base<F>, VR, STAR> &wImag, DistMatrix<Complex<Base<F>>> &Z, int a, int b, SortType sort = UNSORTED)

Compute the eigenpairs of a skew-Hermitian matrix G with indices in the range \(a,a+1,...,b\).

Range-based subset computation

void SkewHermitianEig(UpperOrLower uplo, Matrix<F> &G, Matrix<Base<F>> &wImag, Base<F> a, Base<F> b, SortType sort = UNSORTED)
void SkewHermitianEig(UpperOrLower uplo, DistMatrix<F> &G, DistMatrix<Base<F>, VR, STAR> &wImag, Base<F> a, Base<F> b, SortType sort = UNSORTED)

Compute the eigenvalues of a skew-Hermitian matrix G lying in the half-open interval \((a,b]i\).

void SkewHermitianEig(UpperOrLower uplo, Matrix<F> &G, Matrix<Base<F>> &wImag, Matrix<F> &Z, Base<F> a, Base<F> b, SortType sort = UNSORTED)
void SkewHermitianEig(UpperOrLower uplo, DistMatrix<F> &G, DistMatrix<Base<F>, VR, STAR> &wImag, DistMatrix<F> &Z, Base<F> a, Base<F> b, SortType sort = UNSORTED)

Compute the eigenpairs of a skew-Hermitian matrix G with eigenvalues lying in the half-open interval \((a,b]i\).

Hermitian generalized-definite eigensolvers

The following Hermitian generalized-definite eigenvalue problems frequently appear, where both \(A\) and \(B\) are Hermitian, and \(B\) is additionally positive-definite.

enum HermitianGenDefiniteEigType
enumerator ABX
\[A B x = \lambda x\]
enumerator BAX
\[B A x = \lambda x\]
enumerator AXBX
\[A x = B x \lambda\]

Full spectrum computation

void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, Matrix<F> &A, Matrix<F> &B, Matrix<Base<F>> &w, SortType sort = UNSORTED)
void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<F> &B, DistMatrix<Base<F>, VR, STAR> &w, SortType sort = UNSORTED)

Compute the full set of eigenvalues of a generalized EVP involving the Hermitian matrices A and B, where B is also positive-definite.

void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, Matrix<F> &A, Matrix<F> &B, Matrix<Base<F>> &w, Matrix<Base<F>> &Z, SortType sort = UNSORTED)
void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<F> &B, DistMatrix<Base<F>, VR, STAR> &w, DistMatrix<Base<F>> &Z, SortType sort = UNSORTED)

Compute the full set of eigenpairs of a generalized EVP involving the Hermitian matrices A and B, where B is also positive-definite.

Index-based subset computation

void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, Matrix<F> &A, Matrix<F> &B, Matrix<Base<F>> &w, int a, int b, SortType sort = UNSORTED)
void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<F> &B, DistMatrix<Base<F>, VR, STAR> &w, int a, int b, SortType sort = UNSORTED)

Compute the eigenvalues with indices in the range \(a,a+1,...,b\) of a generalized EVP involving the Hermitian matrices A and B, where B is also positive-definite.

void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, Matrix<F> &A, Matrix<F> &B, Matrix<Base<F>> &w, Matrix<F> &Z, int a, int b, SortType sort = UNSORTED)
void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<F> &B, DistMatrix<Base<F>, VR, STAR> &w, DistMatrix<F> &Z, int a, int b, SortType sort = UNSORTED)

Compute the eigenpairs with indices in the range \(a,a+1,...,b\) of a generalized EVP involving the Hermitian matrices A and B, where B is also positive-definite.

Range-based subset computation

void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, Matrix<F> &A, Matrix<F> &B, Matrix<Base<F>> &w, Base<F> a, Base<F> b, SortType sort = UNSORTED)
void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<F> &B, DistMatrix<Base<F>, VR, STAR> &w, Base<F> a, Base<F> b, SortType sort = UNSORTED)

Compute the eigenvalues lying in the half-open interval \((a,b]\) of a generalized EVP involving the Hermitian matrices A and B, where B is also positive-definite.

void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, Matrix<F> &A, Matrix<F> &B, Matrix<Base<F>> &w, Matrix<F> &Z, Base<F> a, Base<F> b, SortType sort = UNSORTED)
void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<F> &B, DistMatrix<Base<F>, VR, STAR> &w, DistMatrix<F> &Z, Base<F> a, Base<F> b, SortType sort = UNSORTED)

Compute the eigenpairs whose eigenvalues lie in the half-open interval \((a,b]\) of a generalized EVP involving the Hermitian matrices A and B, where B is also positive-definite.

Unitary eigensolver

Not yet written, will likely be based on Ming Gu’s unitary Divide and Conquer algorithm for unitary Hessenberg matrices. The spectral divide and conquer technique described below should suffice in the meantime.

Normal eigensolver

Not yet written, will likely be based on Angelika Bunse-Gerstner et al.’s Jacobi-like method for simultaneous diagonalization of the commuting Hermitian and skew-Hermitian portions of the matrix. The spectral divide and conquer scheme described below should suffice in the meantime.

Schur decomposition

Only a prototype spectral divide and conquer implementation is currently available, though Elemental will eventually also include an implementation of Granat et al.’s parallel QR algorithm.

void Schur(Matrix<F> &A)
void Schur(Matrix<F> &A, Matrix<F> &Q)

Currently defaults to the sequential Hessenberg QR algorithm.

void Schur(DistMatrix<F> &A)
void Schur(DistMatrix<F> &A, DistMatrix<F> &Q)

Currently defaults to the prototype spectral divide and conquer approach.

Hessenberg QR algorithm

void schur::QR(Matrix<F> &A, Matrix<Complex<Base<F>>> &w)
void schur::QR(Matrix<F> &A, Matrix<Complex<Base<F>>> &w, Matrix<F> &Q)

Use a sequential QR algorithm to compute the Schur decomposition.

Spectral divide and conquer

The primary reference for this approach is Demmel et al.’s Fast linear algebra is stable. While the current implementation needs a large number of algorithmic improvements, especially with respect to choosing the Mobius transformations, it tends to succeed on random matrices.

void schur::SDC(Matrix<F> &A, Matrix<Complex<Base<F>>> &w, bool formATR = false, int cutoff = 256, int maxInnerIts = 1, int maxOuterIts = 10, Base<F> relTol = 0)
void schur::SDC(DistMatrix<F> &A, DistMatrix<Complex<Base<F>>, VR, STAR> &w, bool formATR = false, int cutoff = 256, int maxInnerIts = 1, int maxOuterIts = 10, Base<F> relTol = 0)

Compute the eigenvalues of the matrix \(A\) via a spectral divide and conquer process. On exit, the eigenvalues of \(A\) will be stored on its diagonal, and, if formATR was set to true, the upper triangle of \(A\) will be its corresponding upper-triangular Schur factor.

The cutoff controls when the problem is sufficiently small to switch to a sequential Hessenberg QR algorithm, the number of inner iterations is how many attempts to make with the same randomized URV decomposition, and the number of outer iterations is how many random Mobius transformations to try for each spectral split before giving up.

void schur::SDC(Matrix<F> &A, Matrix<Complex<Base<F>>> &w, Matrix<F> &Q, bool formATR = true, int cutoff = 256, int maxInnerIts = 1, int maxOuterIts = 10, Base<F> relTol = 0)
void schur::SDC(DistMatrix<F> &A, DistMatrix<Complex<Base<F>>, VR, STAR> &w, DistMatrix<F> &Q, bool formATR = true, int cutoff = 256, int maxInnerIts = 1, int maxOuterIts = 10, Base<F> relTol = 0)

Attempt to also compute the Schur vectors.

Hermitian SVD

Given an eigenvalue decomposition of a Hermitian matrix \(A\), say

\[A = V \Lambda V^H,\]

where \(V\) is unitary and \(\Lambda\) is diagonal and real. Then an SVD of \(A\) can easily be computed as

\[A = U |\Lambda| V^H,\]

where the columns of \(U\) equal the columns of \(V\), modulo sign flips introduced by negative eigenvalues.

void HermitianSVD(UpperOrLower uplo, Matrix<F> &A, Matrix<Base<F>> &s, Matrix<F> &U, Matrix<F> &V)
void HermitianSVD(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &s, DistMatrix<F> &U, DistMatrix<F> &V)

Return a vector of singular values, \(s\), and the left and right singular vector matrices, \(U\) and \(V\), such that \(A=U \mathrm{diag}(s) V^H\).

void HermitianSVD(UpperOrLower uplo, Matrix<F> &A, Matrix<Base<F>> &s)
void HermitianSVD(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &s)

Return the singular values of \(A\) in s. Note that the appropriate triangle of A is overwritten during computation.

Polar decomposition

Every matrix \(A\) can be written as \(A=QP\), where \(Q\) is unitary and \(P\) is Hermitian and positive semi-definite. This is known as the polar decomposition of \(A\) and can be constructed as \(Q := U V^H\) and \(P := V \Sigma V^H\), where \(A = U \Sigma V^H\) is the SVD of \(A\). Alternatively, it can be computed through (a dynamically-weighted) Halley iteration.

void Polar(Matrix<F> &A)
void Polar(DistMatrix<F> &A)
void Polar(Matrix<F> &A, Matrix<F> &P)
void Polar(DistMatrix<F> &A, DistMatrix<F> &P)

Compute the polar decomposition of \(A\), \(A=QP\), returning \(Q\) within A and \(P\) within P. The current implementation first computes the SVD.

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

Compute the polar decomposition through a Hermitian EVD. Since this is equivalent to a Hermitian sign decomposition (if \(\text{sgn}(0)\) is set to 1), these routines are equivalent to HermitianSign.

polar namespace

int polar::QDWH(Matrix<F> &A, bool colPiv = false, int maxits = 20)
int polar::QDWH(DistMatrix<F> &A, bool colPiv = false, int maxIts = 20)
int hermitian_polar::QDWH(UpperOrLower uplo, Matrix<F> &A, bool colPiv = false, int maxits = 20)
int hermitian_polar::QDWH(UpperOrLower uplo, DistMatrix<F> &A, bool colPiv = false, int maxIts = 20)

Overwrites \(A\) with the \(Q\) from the polar decomposition using a QR-based dynamically weighted Halley iteration. The number of iterations used is returned upon completion. TODO: reference to Yuji’s paper

int polar::QDWH(Matrix<F> &A, Matrix<F> &P, bool colPiv = false, int maxits = 20)
int polar::QDWH(DistMatrix<F> &A, DistMatrix<F> &P, bool colPiv = false, int maxIts = 20)
int hermitian_polar::QDWH(UpperOrLower uplo, Matrix<F> &A, Matrix<F> &P, bool colPiv = false, int maxits = 20)
int hermitian_polar::QDWH(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<F> &P, bool colPiv = false, int maxIts = 20)

Return the full polar decomposition, where \(P\) is HPD.

SVD

Given a general matrix \(A\), the Singular Value Decomposition is the triplet \((U,\Sigma,V)\) such that

\[A = U \Sigma V^H,\]

where \(U\) and \(V\) are unitary, and \(\Sigma\) is diagonal with non-negative entries.

void SVD(Matrix<F> &A, Matrix<Base<F>> &s, Matrix<F> &V)
void SVD(DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &s, DistMatrix<F> &V)

Overwrites A with \(U\), s with the diagonal entries of \(\Sigma\), and V with \(V\).

void SVD(Matrix<F> &A, Matrix<Base<F>> &s)
void SVD(DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &s)

Forms the singular values of \(A\) in s. Note that A is overwritten in order to compute the singular values.

svd namespace

void svd::QRSVD(Matrix<F> &A, Matrix<Base<F>> &s, Matrix<F> &V)

SVD which uses bidiagonal QR algorithm.

void svd::DivideAndConquerSVD(Matrix<F> &A, Matrix<Base<F>> &s, Matrix<F> &V)

SVD which uses a bidiagonal divide-and-conquer algorithm.

void svd::Chan(DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &s, double heightRatio = 1.2)
void svd::Chan(DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &s, DistMatrix<F> &V, double heightRatio = 1.5)

SVD which preprocesses with an initial QR decomposition if the matrix is sufficiently tall relative to its width.

void svd::GolubReinschUpper(DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &s)
void svd::GolubReinschUpper(DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &s, DistMatrix<F> &V)

Computes the singular values (and vectors) of a matrix which is taller than it is wide using the Golub-Reinsch algorithm, though DQDS is used when only the singular values are sought.

void svd::Thresholded(Matrix<F> &A, Matrix<Base<F>> &s, Matrix<F> &V, Base<F> tol = 0, bool relative = false)
void svd::Thresholded(DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &s, DistMatrix<F> &V, Base<F> tol = 0, bool relative = false)

Computes the singular triplets whose singular values are larger than a specified tolerance using the cross-product algorithm. This is often advantageous because tridiagonal eigensolvers tend to enjoy better parallel implementations than bidiagonal SVD’s.