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 .
