The SDC algorithm with inverse free iteration   Next: Spectral Transformation Techniques Up: Parallel Spectral Divide Previous: The SDC algorithm

## The SDC algorithm with inverse free iteration

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  and modified by Bai, Demmel and Gu , 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 .

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 , 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 .   Next: Spectral Transformation Techniques Up: Parallel Spectral Divide Previous: The SDC algorithm

Jack Dongarra
Mon Jan 9 11:07:50 EST 1995