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.
