The eigendecomposition of
an n-by-n real symmetric matrix is the
factorization (
in the complex Hermitian
case), where Z is orthogonal (unitary) and
is real and diagonal,
with
.
The
are the eigenvalues of A, and the columns
of
Z are the eigenvectors . This is also often written
. The eigendecomposition of a symmetric matrix is
computed
by the driver routines PxSYEV and PxSYEVX.
The complex counterparts of these routines, which compute the
eigendecomposition of complex Hermitian matrices, are
the driver routines PxHEEV and PxHEEVX
(see section 3.2.3).
The approximate error
bounds
for the computed eigenvalues
are
The approximate error bounds for the computed eigenvectors
, which bound the acute angles
between the computed eigenvectors and true
eigenvectors
, are
These bounds can be computed by the following code fragment:
EPSMCH = PSLAMCH( ICTXT, 'E' ) * Compute eigenvalues and eigenvectors of A * The eigenvalues are returned in W * The eigenvector matrix Z overwrites A CALL PSSYEV( 'V', UPLO, N, A, IA, JA, DESCA, W, Z, IZ, JZ, $ DESCZ, WORK, LWORK, INFO ) IF( INFO.GT.0 ) THEN PRINT *,'PSSYEV did not converge' ELSE IF( N.GT.0 ) THEN * Compute the norm of A ANORM = MAX( ABS( W( 1 ) ), ABS( W( N ) ) ) EERRBD = EPSMCH * ANORM * Compute reciprocal condition numbers for eigenvectors CALL SDISNA( 'Eigenvectors', N, N, W, RCONDZ, INFO ) DO 10 I = 1, N ZERRBD( I ) = EPSMCH * ( ANORM / RCONDZ( I ) ) 10 CONTINUE END IF
For example, if and
then the eigenvalues, approximate error bounds, and true errors are
as follows.
Further Details: Error Bounds for the Symmetric Eigenproblem
The usual error analysis of the symmetric eigenproblem using ScaLAPACK driver PSSYEV (see subsection 3.2.3) is as follows [101]:
The computed eigendecompositionis nearly the exact eigendecomposition of A+E, namely,
is a true eigendecomposition so that
is orthogonal, where
and
. Here p(n) is a modestly growing function of n. We take p(n)=1 in the above code fragment. Each computed eigenvalue
differs from a true
by at most
Thus, large eigenvalues (those near) are computed to high relative accuracy, while small ones may not be.
The angular difference between the computed unit eigenvector
and a true unit eigenvector
satisfies the approximate bound
ifis small enough. Here,
is the absolute gap between
and the nearest other eigenvalue. Thus, if
is close to other eigenvalues, its corresponding eigenvector
may be inaccurate. The gaps may be easily computed from the array of computed eigenvalues by using subroutine SDISNA . The gaps computed by SDISNA are ensured not to be so small as to cause overflow when used as divisors.
Let
be the invariant subspace spanned by a collection of eigenvectors
, where
is a subset of the integers from 1 to n. Let
be the corresponding true subspace. Then
where
is the absolute gap between the eigenvalues inand the nearest other eigenvalue. Thus, a cluster of close eigenvalues that is far away from any other eigenvalue may have a well-determined invariant subspace
even if its individual eigenvectors are ill-conditioned.
A small possibility exists that PSSYEV will fail to achieve the above error bounds on a heterogeneous network of processors for reasons discussed in section 6.2. On a homogeneous network, PSSYEV is as robust as the corresponding LAPACK routine SSYEV. A future release will attempt to detect heterogeneity and warn the user to use an alternative algorithm.
In contrast to LAPACK, where the same error analysis applied to the
simple and expert drivers, the expert driver PSSYEVX
satisfies slightly weaker error bounds than PSSYEV. The bounds
and
continue to hold,
but the computed eigenvectors
are no longer guaranteed to
be orthogonal to one another. The corresponding LAPACK routine
SSYEVX tries to guarantee orthogonality by reorthogonalizing
computed eigenvectors against one another provided their corresponding
computed eigenvalues are close enough together:
, where the threshold
. If m eigenvalues lie in a cluster satisfying
this closeness criterion, the SSYEVX requires
serial time to
execute. When m is a large fraction of n, this serial bottleneck is
expensive and does not always improve orthogonality.
ScaLAPACK addresses this problem in two ways. First, it lets the user use more or less time and space to perform reorthogonalization, rather than have a fixed criterion. In particular, the user can set the threshold ORTOL used above, decreasing it to make reorthogonalization less frequent or increasing it to reorthogonalize more. Furthermore, since each processor computes a subset of the eigenvectors, ScaLAPACK permits reorthogonalization only with the local eigenvectors; that is, no communication is allowed. Hence, if a cluster of eigenvalues is small enough for the corresponding eigenvectors to fit on one processor, the same reorthogonalization will be done as in LAPACK. The user can supply more or less workspace to limit the size of a cluster on one processor. Hence, at one extreme, with a large cluster and lots of workspace, the algorithm will be essentially equivalent to SSYEVX. At the other extreme, with all small clusters or little workspace supplied, the algorithm will be perfectly load balanced and perform minimal communication to compute the eigenvectors.
The second way ScaLAPACK will deal with reorthogonalization is to introduce a new algorithm [103, 102, 44] that requires nearly no reorthogonalization to compute orthogonal eigenvectors in a fully parallel way. This algorithm will be introduced in future ScaLAPACK and LAPACK releases.
In the special case of a real symmetric tridiagonal matrix T, the eigenvalues can sometimes be computed much more accurately. PxSYEV (and the other symmetric eigenproblem drivers) computes the eigenvalues and eigenvectors of a dense symmetric matrix by first reducing it to tridiagonal form T and then finding the eigenvalues and eigenvectors of T. Reduction of a dense matrix to tridiagonal form T can introduce additional errors, so the following bounds for the tridiagonal case do not apply to the dense case.
The eigenvalues of the symmetric tridiagonal matrix T may be computed with small componentwise relative backward error () by using subroutine PxSTEBZ (section 3.3.4). To compute tighter error bounds for the computed eigenvalues
we must make some assumptions about T. The bounds discussed here are from [15]. Suppose T is positive definite, and write T=DHD where
and
. Then the computed eigenvalues
can differ from true eigenvalues
by
where p(n) is a modestly growing function of n. Thus, ifis moderate, each eigenvalue will be computed to high relative accuracy, no matter how tiny it is.