next up previous contents index
Next: Error Bounds for the Up: Accuracy and Stability Previous: Error Bounds for Linear

Error Bounds for the Symmetric Eigenproblem


The eigendecomposition  of an n-by-n real symmetric matrix is the factorization tex2html_wrap_inline13904 (tex2html_wrap_inline18618 in the complex Hermitian case), where Z is orthogonal (unitary) and tex2html_wrap_inline18622 is real and diagonal, with tex2html_wrap_inline18624. The tex2html_wrap_inline18626 are the eigenvalues  of A, and the columns tex2html_wrap_inline18630 of Z are the eigenvectors . This is also often written tex2html_wrap_inline18634. 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 tex2html_wrap_inline18636 are
The approximate error bounds for the computed eigenvectors tex2html_wrap_inline18638, which bound the acute angles between the computed eigenvectors and true eigenvectors tex2html_wrap_inline18630, are    
These bounds can be computed by the following code fragment:

*     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 ) ) )
*        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 tex2html_wrap_inline18242 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 eigendecomposition tex2html_wrap_inline18680 is nearly the exact eigendecomposition of A+E, namely, tex2html_wrap_inline18684 is a true eigendecomposition so that tex2html_wrap_inline18686 is orthogonal, where tex2html_wrap_inline18688 and tex2html_wrap_inline18690. Here p(n) is a modestly growing function of n. We take p(n)=1 in the above code fragment. Each computed eigenvalue tex2html_wrap_inline18698 differs from a true tex2html_wrap_inline18626 by at most
Thus, large eigenvalues (those near tex2html_wrap_inline18702) are computed to high relative accuracy,   while small ones may not be.  

The angular difference between the computed unit eigenvector tex2html_wrap_inline18704 and a true unit eigenvector tex2html_wrap_inline18630 satisfies the approximate bound  
if tex2html_wrap_inline18543 is small enough. Here, tex2html_wrap_inline18712 is the absolute gap   between tex2html_wrap_inline18626 and the nearest other eigenvalue. Thus, if tex2html_wrap_inline18626 is close to other eigenvalues, its corresponding eigenvector tex2html_wrap_inline18630 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 tex2html_wrap_inline17614 be the invariant subspace spanned by a collection of eigenvectors tex2html_wrap_inline18722, where tex2html_wrap_inline18724 is a subset of the integers from 1 to n. Let tex2html_wrap_inline17602 be the corresponding true subspace. Then
is the absolute gap between the eigenvalues in tex2html_wrap_inline18724 and 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 tex2html_wrap_inline17614 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 tex2html_wrap_inline18736 and tex2html_wrap_inline18738 continue to hold, but the computed eigenvectors tex2html_wrap_inline18638 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: tex2html_wrap_inline18742, where the threshold tex2html_wrap_inline18744. If m eigenvalues lie in a cluster satisfying this closeness criterion, the SSYEVX requires tex2html_wrap_inline18748 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     (tex2html_wrap_inline18764) by using subroutine PxSTEBZ (section 3.3.4).    To compute tighter error bounds for the computed eigenvalues tex2html_wrap_inline18698 we must make some assumptions about T. The bounds discussed here are from [15]. Suppose T is positive definite, and write T=DHD where tex2html_wrap_inline18774 and tex2html_wrap_inline18776. Then the computed eigenvalues tex2html_wrap_inline18698 can differ from true eigenvalues tex2html_wrap_inline18626 by
where p(n) is a modestly growing function of n. Thus, if tex2html_wrap_inline18786 is moderate, each eigenvalue will be computed to high relative accuracy,   no matter how tiny it is.

next up previous contents index
Next: Error Bounds for the Up: Accuracy and Stability Previous: Error Bounds for Linear

Susan Blackford
Tue May 13 09:21:01 EDT 1997