The above algorithm needs an explicit matrix inverse. This could cause numerical instability when the matrix is ill-conditioned. The following algorithm, originally due to Godunov, Bulgakov and Malyshev [44][10][30] and modified by Bai, Demmel and Gu [7], eliminates the need for the matrix inverse, and divides the spectrum along the unit circle instead of the imaginary axis. We first describe the algorithm, and then briefly explain why it works.

ALGORITHM 2 (THE SDC ALGORITHM WITH INVERSE FREE ITERATION)

Let , ;

For until convergence or do

; ~ (QR decomposition)

if ), , exit

End for;

Compute ; ~ (rank revealing QR decomposition)

Let ;

Compute ;

Compute .

As in Algorithm 1, we need to choose a stopping criterion in the inner loop, as well as a limit on the maximum number of iterations. On convergence, the eigenvalues of are the eigenvalues of inside the unit disk , and the eigenvalues of are the eigenvalues of outside . It is assumed that no eigenvalues of are on the unit circle. As with Algorithm 1, the quantity measures the backward stability.

To illustrate how the algorithm works we will assume that all matrices we want to invert are invertible. From the inner loop of the algorithm, we see that

so or . Therefore

so the algorithm is simply repeatedly squaring the eigenvalues, driving the ones inside the unit circle to 0 and those outside to . Repeated squaring yields quadratic convergence. This is analogous to the sign function iteration where computing is equivalent to taking the Cayley transform of , squaring, and taking the inverse Cayley transform. Further explanation of how the algorithm works can be found in [7].

An attraction of this algorithm is that it can equally well deal with the generalized nonsymmetric eigenproblem , provided the problem is regular, i.e. is not identically zero. One simply has to start the algorithm with instead of .

Regarding the QR decomposition in the inner loop, there is no need to form the entire unitary matrix in order to get the submatrices and . Instead, we can compute the QR decomposition of the matrix , which leaves stored implicitly as Householder vectors in the lower triangular part of the matrix and another dimensional array. We can then apply - without computing it - to the matrix to obtain the desired matrices and .

We now show how to compute in the rank revealing QR decomposition
of without computing
the explicit inverse and subsequent products.
This will yield the ultimate *inverse free* algorithm.
Recall that for our purposes, we only need the unitary factor and the
rank of .
It turns out that by using the generalized QR decomposition technique
developed in [2][45],
we can get the desired information without computing
. In fact, in order to compute the QR decomposition with
pivoting of , we first compute the QR decomposition
with pivoting of the matrix :

and then we compute the RQ factorization of the matrix :

From (2.7) and (2.8), we have . The is the desired unitary factor. The rank of is also the rank of the matrix .

Mon Jan 9 11:07:50 EST 1995