Matrix decompositions¶
Hermitian eigensolver¶
Elemental provides a collection of routines for both full and partial solution of the Hermitian eigenvalue problem
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:
Compute all eigenvalues or eigenpairs, \(\{\lambda_i\}_{i=0}^{n-1}\) or \(\{(z_i,\lambda_i)\}_{i=0}^{n-1}\).
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\).
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, 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
This fact implies a fast method for solving skew-Hermitian eigenvalue problems:
Form \(iG\) in \(O(n^2)\) work (switching to complex arithmetic in the real case)
Run a Hermitian eigensolve on \(iG\), yielding \(iG=Z \Lambda Z^H\).
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\]
-
enumerator
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, 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¶
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
where \(V\) is unitary and \(\Lambda\) is diagonal and real. Then an SVD of \(A\) can easily be computed as
where the columns of \(U\) equal the columns of \(V\), modulo sign flips introduced by negative eigenvalues.
-
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
(DistMatrix<F> &A)¶
-
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
(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
(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
where \(U\) and \(V\) are unitary, and \(\Sigma\) is diagonal with non-negative entries.
-
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
(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.