Numerical Stability

Robust implementations of a practical QR algorithm, such as are available in LAPACK [12], compute an approximate Schur form for a matrix $A$ that satisfies

$$(A+E)\hat{V} = \hat{V} \hat{R},$$ (128)

where $\hat{R}$ is upper triangular, $\Vert \hat{V}^\ast \hat{V} - I\Vert \approx \epsilon_M$, and $\Vert E\Vert \approx \epsilon_M\Vert A\Vert$, where $\epsilon_M$ is machine precision (unit roundoff). The fact that $\Vert E\Vert \approx \epsilon_M\Vert A\Vert$ is a direct result of the exclusive use of unitary transformations during the iteration. This is what makes the QR algorithm backward stable. In other words, the QR method computes the Schur form of a matrix $A+E$ that is near $A$.

Backward stability is a very reassuring quality for any numerical algorithm, but it is important to keep in mind that backward stability alone does not imply accurate answers. The accuracy of the computed eigenvalues and eigenvectors depends upon the sensitivity (conditioning) of the eigensystem of $A$. A backward stable algorithm produces accurate answers for well-conditioned problems. The reader is referred to §7.13 for a more detailed discussion of this issue.

As explained in the last section, the IRAM is a truncated QR iteration. In most implementations of the QR method, Householder transformations are used to obtain the initial reduction to Hessenberg form. The Arnoldi procedure, if extended to $n$ steps, is simply another way to obtain this reduction. The introduction of the orthogonal correction [96] allows the Arnoldi procedure to produce a reduction of the same numerical quality as the Householder reduction but it is more suitable to the large scale setting. With this correction, the computed Arnoldi vectors are unitary to within machine precision $(\Vert V_k^*V_k - I_k \Vert \approx \epsilon_M)$, and since the restarting mechanism involves the same unitary transformations used in the QR mechanism, the IRAM is also backward stable.

At any stage of the IRAM iteration, we have

$$(A + E_o) V_k = V_k H_k + f_k e_k^*,$$

with $\Vert E_o\Vert \approx \Vert A\Vert\epsilon_M$. If we are successful in making $\Vert f_k\Vert < \Vert A\Vert\epsilon_M$ then

$$(A + E) V_k = V_k H_k,$$

with $E = E_o + f_k (V_k e_k)^*$. If $H_k Z_k = Z_k R_k$ is a Schur decomposition of $H_k$ then

$$(A + E) \hat{V}_k = \hat{V}_k R_k \ \ \mbox{with} \ \ \hat{V}_k = V_k Z_k$$

is an exact partial Schur decomposition of a nearby matrix $A+E$.

It is more likely that convergence will be realized with a small last component of an eigenvector of $H_m$ ($m = k+p$) at some stage. We still have

$$(A + E_1) V_m = V_m H_m + f_m e_m^*,$$

with $E_1$ small relative to $\Vert A\Vert$. If $H_m y = y \theta$ where $\eta = e_m^*y$ is small ( $\Vert f_m\Vert \vert\eta\vert < \epsilon_U \Vert H_m \Vert$, say), then (A + E) x = x , with $x = V_m y$ and $E = E_1 + f_m (e_m^*y) x^*$ (assuming $\Vert y\Vert = 1$). Thus, in any case, converged approximate eigenpairs $(x, \theta)$ are exact for a nearby matrix $A+E$.

Through the use of the deflation techniques that we shall present here in §7.6.6, it is possible to use unitary transformations to obtain a backward stability statement with respect to a partial Schur decomposition for a set of $k$ converged eigenvalues. If there are $k$ eigenvalues of $H_m$ that have converged, it is possible to construct a unitary matrix $Q$ such that the leading $k \times k$ submatrix $R_k$ of $Q^* H_m Q$ is upper triangular and

$$( A + \hat{E}) \hat{V}_k = \hat{V}_k R_k , \ \ \mbox{with} \ \ \hat{V}_k = V_m Q,$$

where $\hat{E}$ is small relative to $\Vert A\Vert$. A discussion of this is given in §7.6.6, and more detail is available in [420].