The first SDC algorithm uses the matrix sign function, which was introduced by Roberts [46] for solving the algebraic Riccati equation. However, it was soon extended to solving the spectral decomposition problem [8]. More recent studies may be found in [6][42][11].

The matrix sign function, , of a matrix with no eigenvalues on the imaginary axis can be defined via the Jordan canonical form of (2.1), where the eigenvalues of are in the open right half plane , and the eigenvalues of are in the open left half plane . Then is

It is easy to see that the matrix

is the spectral projector onto the invariant subspace corresponding to the eigenvalues of in . is the number of the eigenvalues of in . is the spectral projector corresponding to the eigenvalues of in .

Now let be the rank revealing QR decomposition of the projector . Then yields the desired spectral decomposition (2.3), where the eigenvalues of are the eigenvalues of in , and the eigenvalues of are the eigenvalues of in .

Since the matrix sign function, , satisfies the matrix equation , we can use Newton's method to solve this matrix equation and obtain the following simple iteration:

The iteration is globally and ultimately quadratically convergent with , provided has no pure imaginary eigenvalues [35][46]. The iteration fails otherwise, and in finite precision, the iteration could converge slowly or not at all if is ``close'' to having pure imaginary eigenvalues.

There are many ways to improve the accuracy and convergence rate of this basic iteration [37][33][12]. For example, if , we may use the so called Newton-Schulz iteration

to avoid the use of the matrix inverse. Although it requires twice as many flops, it is more efficient whenever matrix multiply is at least twice as efficient as matrix inversion. The Newton-Schulz iteration is also quadratically convergent provided that . A hybrid iteration might begin with Newton iteration until and then switch to Newton-Schulz iteration (we discuss the performance of one such hybrid later).

Hence, we have the following algorithm which divides the spectrum along the pure imaginary axis.

ALGORITHM 1 (THE SDC ALGORITHM WITH NEWTON ITERATION)

Let ;

For until convergence or do

;

if , , exit

End for;

Compute ; ~ (rank revealing QR decomposition)

Let ;

Compute ;

Compute .

Here is the stopping criterion for the Newton iteration (say, , where is the machine precision), and limits the maximum number of iterations (say ). On return, the generally nonzero quantity measures the backward stability of the computed decomposition, since by setting to zero and so decoupling the problem into and , a backward error of is introduced.

For simplicity, we just use the QR decomposition with column pivoting to reveal rank, although more sophisticated rank-revealing schemes exist [51][34][31][14].

All the variations of the Newton iteration with global convergence still need to compute the inverse of a matrix explicitly in one form or another. Dealing with ill-conditioned matrices and instability in the Newton iteration for computing the matrix sign function and the subsequent spectral decomposition is discussed in [4][6][11] and the references therein.

Mon Jan 9 11:07:50 EST 1995