Further Details: Error Bounds for the Symmetric Eigenproblem

The usual error analysis of the symmetric eigenproblem (using any LAPACK routine in subsection 2.3.4 or any EISPACK routine) is as follows [85]:

The computed eigendecomposition $\hat{Z} \hat{\Lambda} \hat{Z}^T$ is nearly the exact eigendecomposition of A+E, i.e., $A+E = (\hat{Z}+ \delta \hat{Z}) \hat{\Lambda} (\hat{Z}+ \delta \hat{Z})^T$ is a true eigendecomposition so that $\hat{Z}+ \delta \hat{Z}$ is orthogonal, where $\Vert E\Vert _2 / \Vert A\Vert _2 \leq p(n) \epsilon$ and $\Vert \delta \hat{Z} \Vert _2 \leq p(n) \epsilon$. Here p(n) is a modestly growing function of n. We take p(n)=1 in the above code fragment. Each computed eigenvalue $\hat{\lambda}_i$ differs from a true $\lambda_i$ by at most

$$\vert \hat{\lambda}_i - \lambda_i \vert \leq p(n) \cdot \epsilon \cdot \Vert A\Vert _2 = {\tt EERRBD} \; .$$

Thus large eigenvalues (those near $\max_i \vert \lambda_i \vert = \Vert A\Vert _2$) are computed to high relative accuracy and small ones may not be.

There are two questions to ask about the computed eigenvectors: ``Are they orthogonal?'' and ``How much do they differ from the true eigenvectors?'' The answer to the first question is that except for eigenvectors computed by computational routine xSTEIN (which is called by drivers with names ending in -EVX when only a subset of the eigenvalues and eigenvectors are requested), the computed eigenvectors are always nearly orthogonal to working precision, independent of how much they differ from the true eigenvectors. In other words

$$\vert\hat{z}_i^T \hat{z}_j \vert = O( \epsilon )$$

for $i \neq j$. xSTEIN almost always returns orthogonal eigenvectors, but can occasionally fail when eigenvalues are tightly clustered.

Here is the answer to the second question about eigenvectors. The angular difference between the computed unit eigenvector $\hat{z}_i$ and a true unit eigenvector z_i satisfies the approximate bound

$$\theta ( \hat{z}_i , z_i ) \mathrel{\raisebox{-.75ex}{$\math... ...{p(n) \epsilon \Vert A\Vert _2}{{\rm gap}_i} = {\tt ZERRBD}(i)$$

if $p(n) \epsilon$ is small enough. Here ${\rm gap}_i = \min_{j \neq i} \vert \lambda_i - \lambda_j \vert$ is the absolute gap between $\lambda_i$ and the nearest other eigenvalue. Thus, if $\lambda_i$ is close to other eigenvalues, its corresponding eigenvector z_i may be inaccurate. The gaps may be easily computed from the array of computed eigenvalues using subroutine SDISNA. The gaps computed by SDISNA are ensured not to be so small as to cause overflow when used as divisors.

Let $\hat{\cal S}$ be the invariant subspace spanned by a collection of eigenvectors $\{\hat{z}_i \, , \, i \in {\cal I}\}$, where $\cal I$ is a subset of the integers from 1 to n. Let $\cal S$ be the corresponding true subspace. Then

$$\theta ( {\hat{\cal S}}, {\cal S}) \mathrel{\raisebox{-.75ex... ...<}$}}\frac{p(n) \epsilon \Vert A\Vert _2} {{\rm gap}_{\cal I}}$$

where

$${\rm gap}_{\cal I} = \min_{i \in {\cal I} \atop j \not\in {\cal I}} \vert \lambda_i - \lambda_j \vert$$

is the absolute gap between the eigenvalues in $\cal I$ and the nearest other eigenvalue. Thus, a cluster of close eigenvalues which is far away from any other eigenvalue may have a well determined invariant subspace $\hat{\cal S}$ even if its individual eigenvectors are ill-conditioned^4.3.

In the special case of a real symmetric tridiagonal matrix T, the eigenvalues and eigenvectors can be computed much more accurately. xSYEV (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 T may be computed with small componentwise relative backward error ($O(\epsilon)$) by using subroutine xSTEBZ (subsection 2.4.4) or driver xSTEVX (subsection 2.3.4). If T is also positive definite, they may also be computed at least as accurately by xPTEQR or xSTEGR (subsection 2.4.4). To compute error bounds for the computed eigenvalues $\hat{\lambda}_i$ we must make some assumptions about T. The bounds discussed here are from [14]. Suppose T is positive definite, and write T=DHD where $D = {\mbox {\rm diag}}( t_{11}^{1/2} , \ldots , t_{nn}^{1/2})$ and h_ii= 1. Then the computed eigenvalues $\hat{\lambda}_i$ can differ from true eigenvalues $\lambda_i$ by

$$\vert \hat{\lambda}_i - \lambda_i \vert \leq p(n) \cdot \epsilon \cdot \kappa_2(H) \cdot \lambda_i$$

where p(n) is a modestly growing function of n. Thus if $\kappa_2 (H)$ is moderate, each eigenvalue will be computed to high relative accuracy, no matter how tiny it is. The eigenvectors $\hat{z}_i$ computed by xPTEQR can differ from true eigenvectors z_i by at most about

$$\theta ( \hat{z}_i , z_i ) \mathrel{\raisebox{-.75ex}{$\math... ...\frac{p(n) \cdot \epsilon \cdot \kappa_2 (H)} {{\rm relgap}_i}$$

if $p(n) \epsilon$ is small enough, where ${\rm relgap}_i = \min_{j \neq i} \vert \lambda_i - \lambda_j \vert/ ( \lambda_i + \lambda_j )$ is the relative gap between $\lambda_i$ and the nearest other eigenvalue. Since the relative gap may be much larger than the absolute gap, this error bound may be much smaller than the previous one. Independent of the difference between the computed and true eigenvectors, the computed eigenvectors are orthogonal to nearly full precision, i.e. $\vert \hat{z}_i^T \hat{z}_j \vert = O( \epsilon )$ for $i \neq j$.

$\kappa_2 (H)$ could be computed by applying xPTCON (subsection 2.4.1) to H. The relative gaps are easily computed from the array of computed eigenvalues.