Further Details: Error Bounds for the Generalized Symmetric Definite Eigenproblem

The error analysis of the driver routines for the generalized symmetric definite eigenproblem goes as follows. In all cases ${\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.

1.

$A - \lambda B$. The computed eigenvalues $\hat{\lambda}_i$ can differ from true eigenvalues $\lambda_i$ by at most about

$$\vert\hat{\lambda}_i - \lambda_i \vert \mathrel{\raisebox{-.... ...B) \cdot \vert \hat{\lambda}_i \vert ) = {\tt EERRBD}(i) \; .$$

The angular difference between the computed eigenvector $\hat{z}_i$ and a true eigenvector z_i is

$$\theta ( \hat{z_i} , z_i ) \mathrel{\raisebox{-.75ex}{$\math... ... \hat{\lambda}_i \vert} { {\rm gap}_i } = {\tt ZERRBD}(i) \; .$$

2.

$AB- \lambda I$ or $BA- \lambda I$. The computed eigenvalues $\hat{\lambda}_i$ can differ from true eigenvalues $\lambda_i$ by at most about

$$\vert\hat{\lambda}_i - \lambda_i \vert \mathrel{\raisebox{-.... ...B) \cdot \vert \hat{\lambda}_i \vert ) = {\tt EERRBD}(i) \; .$$

The angular difference between the computed eigenvector $\hat{z}_i$ and a true eigenvector z_i is

$$\theta ( \hat{z_i} , z_i ) \mathrel{\raisebox{-.75ex}{$\math... ... { {\rm gap}_i } + \kappa_2 (B) \right) = {\tt ZERRBD}(i) \; .$$

The code fragments above replace p(n) by 1, and makes sure neither RCONDB nor RCONDZ is so small as to cause overflow when used as divisors in the expressions for error bounds.

These error bounds are large when B is ill-conditioned with respect to inversion ($\kappa_2 (B)$ is large). It is often the case that the eigenvalues and eigenvectors are much better conditioned than indicated here. We mention three ways to get tighter bounds. The first way is effective when the diagonal entries of B differ widely in magnitude^4.1:

1.

$A - \lambda B$. Let $D = {\mbox {\rm diag}}(b_{11}^{-1/2} , \ldots , b_{nn}^{-1/2})$ be a diagonal matrix. Then replace B by DBD and A by DAD in the above bounds.

2.

$AB- \lambda I$ or $BA- \lambda I$. Let $D = {\mbox {\rm diag}}(b_{11}^{-1/2} , \ldots , b_{nn}^{-1/2})$ be a diagonal matrix. Then replace B by DBD and A by D^-1AD^-1 in the above bounds.

The second way to get tighter bounds does not actually supply guaranteed bounds, but its estimates are often better in practice. It is not guaranteed because it assumes the algorithm is backward stable, which is not necessarily true when B is ill-conditioned. It estimates the chordal distance between a true eigenvalue $\lambda_i$ and a computed eigenvalue $\hat{\lambda}_i$:

$$\chi ( \hat{\lambda}_i , \lambda_i ) = \frac{\vert \hat{\lam... ... {\sqrt{1 + \hat{\lambda}_i^2} \cdot \sqrt{1 + \lambda_i^2} } .$$

To interpret this measure we write $\lambda_i = \tan \theta$ and $\hat{\lambda}_i = \tan \hat{\theta}$. Then $\chi ( \hat{\lambda}_i , \lambda_i ) = \vert \sin ( \hat{\theta} - \theta ) \vert$. In other words, if $\hat{\lambda}_i$ represents the one-dimensional subspace $\hat{\cal S}$ consisting of the line through the origin with slope $\hat{\lambda}_i$, and $\lambda_i$ represents the analogous subspace $\cal S$, then $\chi ( \hat{\lambda}_i , \lambda_i )$ is the sine of the acute angle $\theta (\hat{\cal S}, {\cal S})$ between these subspaces. Thus $\chi$ is bounded by one, and is small when both arguments are large^4.2. It applies only to the first problem, $A - \lambda B$:

Suppose a computed eigenvalue $\hat{\lambda}_i$ of $A - \lambda B$ is the exact eigenvalue of a perturbed problem $(A+E) - \lambda (B+F)$. Let x_i be the unit eigenvector (|x_i|₂=1) for the exact eigenvalue $\lambda_i$. Then if |E| is small compared to |A|, and if |F| is small compared to |B|, we have

$$\chi ( \hat{\lambda}_i , \lambda_i ) \mathrel{\raisebox{-.75... ...Vert F\Vert}{\sqrt{(x_i^H Ax_i )^2 + (x_i^H Bx_i )^2}} \; \; .$$

Thus $1/{\sqrt{(x_i^H Ax_i )^2 + (x_i^H Bx_i )^2}}$ is a condition number for eigenvalue $\lambda_i$.

The third way applies only to the first problem $A - \lambda B$, and only when A is positive definite. We use a different algorithm:

1.

Compute the Cholesky factorization of A = U^T_A U_A, using xPOTRF.

2.

Compute the Cholesky factorization of B = U^T_B U_B, using xPOTRF.

3.

Compute the generalized singular value decomposition of the pair U_A, U_B using xTGSJA. The squares of the generalized singular values are the desired eigenvalues.

See sections 2.3.5.3 and  2.4.9 for a discussion of the generalized singular value decomposition, and section 4.12 for a discussion of the relevant error bound. This approach can give a tighter error bound than the above bounds when B is ill conditioned but A+B is well-conditioned.

Other yet more refined algorithms and error bounds are discussed in [14,95,103].