Using the Matrix Sign Function to Compute Invariant Subspaces.
The matrix sign function has several applications in system theory and matrix computations. However, the numerical behavior of the matrix sign function, and its associated divide-and-conquer algorithm for computing invariant subspaces, are still not completely understood. In this paper, we present a new perturbation theory for the matrix sign function,  the conditioning of its computation, the numerical stability of the divide-and-conquer algorithm, and iterative refinement schemes. Numerical examples are also presented. An extension of the matrix-sign-function-based algorithm to compute left and right deflating subspaces for a regular pair of matrices is also described.
Introduction
. Since the matrix sign function was introduced in early 1970s,
it has been the subject of numerous studies and used in many applications. For
example, see [30, 31, 11, 26, 23] and references therein. Our main interest here is
to use the matrix sign function to build parallel algorithms for computing invariant
subspaces of nonsymmetric matrices, as well as their associated eigenvalues. It is
a challenge to design a parallel algorithm for the nonsymmetric eigenproblem that
uses coarse grain parallelism effectively, scales for larger problems on larger machines,
does not waste time dealing with the parts of the spectrum in which the user is not
interested, and deals with highly nonnormal matrices and strongly clustered spectra.
In the work of [2], after reviewing the existing approaches, we proposed a design of a
parallel nonsymmetric eigenroutine toolbox, which includes the basic building blocks
(such as LU factorization, matrix inversion and the matrix sign function), standard
eigensolver routines (such as the QR algorithm) and new algorithms (such spectral
divide-and-conquer using the sign function). We discussed how these tools could be
used in different combinations on different problems and architectures, for extracting
all or some of the eigenvalues of a nonsymmetric matrix, and/or their corresponding
invariant subspaces. Rather than using "black box" eigenroutines such as provided
by EISPACK [32, 21] and LAPACK [1], we expect the toolbox approach to allow us
more flexibility in developing efficient problem-oriented eigenproblem solvers on high
performance machines, especially on parallel distributed memory machines.
However, the numerical accuracy and stability of the the matrix sign function and
divide-and-conquer algorithms based on it are poorly understood. In this paper, we
will address these issues. Much of this work also appears in [3].
Let us first restate some of basic definitions and ideas to establish notation. The
matrix sign function of a matrix A is defined as follows [30]: Let
be the Jordan canonical form of a matrix A 2 C n\Thetan , where the eigenvalues of J+ lie
in the open right half plane (C+ ) and those of J \Gamma lie in the open left half plane (C \Gamma ).
Published at SIAM J. Mat. Anal. Appl., Vol.19, pp.205-225, 1998
y Department of Mathematics, University of Kentucky, Lexington, KY 40506 (bai@ms.uky.edu).
z Computer Science Division and Mathematics Department, University of California, Berkeley,
Z. BAI AND J. DEMMEL
Then the matrix sign function of A is:
We assume that no eigenvalue of A lies on the imaginary axis; otherwise, sign(A)
is not defined. It is easy to show that the spectral projection corresponding to the
eigenvalues of A in the open right and left half planes are P
respectively. Let the leading columns of an orthogonal matrix Q span the range space
of P+ (for example, Q may be computed by the rank-revealing QR decomposition of
the spectral decomposition
(1)
where are the eigenvalues of A in C+ , and -(A 22 ) are the eigenvalues of A
in C \Gamma . The algorithm proceeds in a divide-and-conquer fashion, by computing the
eigenvalues of A 11 and A 22 .
Rather than using the Jordan canonical form to compute sign(A), it can be shown
that sign(A) is the limit of the following Newton iteration
A k+1 =2
(2)
The iteration is globally and ultimately quadratic convergent. There exist different
scaling schemes to speedup the convergence of the iteration, and make it more suitable
for parallel computation. By computing the matrix sign function of a M-obius
transformation of A, the spectrum can be divided along arbitrary lines and circles,
rather than just along the imaginary axis. See the report [2] and the references therein
for more details.
Unfortunately, in finite precision arithmetic, the ill conditioning of a matrix A k
with respect to inversion and rounding errors, may destroy the convergence of the
Newton iteration (2), or cause convergence to the wrong answer. Consequently, the
left bottom corner block of the matrix Q T AQ in (1) may be much larger than
where u denotes machine precision. This means that it is not numerically stable to
approximate the eigenvalues of A by the eigenvalues of A 11 and A 22 , as we would like.
In this paper, we will first study the perturbation theory of the matrix sign func-
tion, its conditioning, and the numerical stability of the overall divide-and-conquer
algorithm based on the matrix-sign function. We realize that it is very difficult to give
a complete and clear analysis. We only have a partial understanding of when we can
expect the Newtion iteration to converge, and how accurate it is. In a coarse analysis,
we can also bound the condition numbers of intermediate matrices in the Newton iter-
ation. Artificial and possibly very pathological test matrices are constructed to verify
our theoretical analysis. Besides these artificial tests, we also test a large number of
eigenvalue problems of random matrices, and a few eigenvalue problems from appli-
cations, such as electrical power system analysis, numerical simulation of chemical
reactions, and areodynamics stability analysis. Through these examples, we conclude
that the most bounds for numerical sensitivity and stability of matrix sign function
computation and its based algorithms are reachable for some very pathological cases,
but they are often very pessimistic. The worst cases happen rarely.
In addition, we discuss iterative refinement of an approximate invariant subspace,
and outline an extension of the matrix sign function based algorithms to compute both
left and right deflating subspaces for a regular matrix pencil A \Gamma -B.
MATRIX SIGN FUNCTION FOR COMPUTING INVARIANT SUBSPACES 3
The rest of this paper is organized as the following: Section 2 presents a new
perturbation bound for the matrix sign function. Section 3 discusses the numerical
conditioning of the matrix sign function. The backward error analysis of computed
invariant subspace and remarks on the matrix sign function based algorithm versus the
QR algorithm are presented in section 4. Section 5 presents some numerical examples
for the analysis of sections 2, 3 and 4. Section 6 describes the iteration refinement
scheme to improve an approximate invariant subspace. Section 7 outlines an extension
of the matrix sign function based algorithms for the generalized eigenvalue problem.
Concluding remarks are presented in section 8.
2. A perturbation bound for the matrix sign function. When a matrix A
has eigenvalues on pure imaginary axis, its matrix sign function is not defined. In other
words, the set of ill-posed problems for the matrix sign function, is the set of matrices
with at least one pure-imaginary eigenvalue. Computationally, we have observed that
when there are the eigenvalues of A close the pure imaginary axis, the Newton iteration
and its variations are very slowly convergent, and may be misconvergent. Moreover,
even when the iteration converges, the error in the computed matrix sign function
could be too large to use. It is desirable to have a perturbation analysis of the matrix
sign function related to the distance from A to the nearest ill-posed problem.
Perturbation theory and condition number estimation of the matrix sign function
are discussed in [25, 23, 29]. However, none of the existing error bounds explicitly
reveals the relationship between the sensitivity of the matrix sign function and the
distance to the nearest ill-posed problem. In this section, we will derive a new perturbation
bound which explicitly reveals such relationship. We will denote all the eigen-values
of A with positive real part by -+ (A), i.e.,
denotes the smallest singular value of A. In addition, we recall the well-known
inequality
is the matrix 2-norm.
Theorem 2.1. Suppose A has no pure imaginary and zero eigenvalues,
is a perturbation of A and ffl j kffiAk. Let
Then
3:
Furthermore, if
then
4 Z. BAI AND J. DEMMEL
O
Re
Im
r
-r
r
Fig. 1. The semi-circle \Gamma
Proof. We only prove the bound (7). The bound (5) can be proved by using a
similar technique. Following Roberts [30] (or Kato [24]), the matrix sign function can
also be defined using Cauchy integral representation:
where
Z
\Gamma is any simple closed curve with positive direction enclosing -+ (A). sign
the spectral projector for -+ (A). Here, without loss of generality, we take \Gamma to be a
semi-circle with radius Figure 1). From the definition
(8) of sign(A), it is seen that to study the stability of the matrix sign function of A
to the perturbation ffi A, it is sufficient to just study the sensitivity of the projection
be the projection corresponding to -+ A), from the
condition (6), no eigenvalues of A are perturbed across or on the pure imaginary axis,
and the semi-circle \Gamma also encloses -+ Therefore we have
Z
Z r
\Gammar
Z -=2
\Gamma-=2
[(re
where the first integral, denoted I 1 , is the integral over the straight line of the semi-circle
\Gamma, the second integral, denoted I 2 , is the integral over the curved part of the
semi-circle \Gamma. Now, by taking the spectral norm of the first integral term, and noting
the definition of !, the condition (6) and the inequality (3), we have
Z r
\Gammar
MATRIX SIGN FUNCTION FOR COMPUTING INVARIANT SUBSPACES 5
Z r
\Gammar
Z r
\Gammar
Z r
\Gammar
By taking the spectral norm of the second integral term I 2 , we have
Z -=2
\Gamma-=2
k(re
Z -=2
\Gamma-=2
re i'
re i'
r
r
where the third inequality follows from (3) and the fourth from the choice of the
radius r of the semi-circle \Gamma. The desired bound (7) follows from the bounds on kI 1 k
and kI 2 k and the identity
A few remarks are in order:
1. In the language of pseudospectra [35], the condition (6) means that the kffiAk-
pseudospectra of A do not cross the pure imaginary axis.
2. From the perturbation bound (7), we see that the stability of the matrix sign
function to the perturbation requires not only the kffiAk-pseudospectra of the
A to be bounded away from the pure imaginary axis, but also
A to
be small (recall that dA is the distance from A to the nearest matrix with a
pure-imaginary eigenvalue).
3. It is natural to take
A as the condition number of the matrix sign
function. Algorithms for computing dA and related problems can be found
in [14, 9, 8, 12].
4. The bound (7) is similar to the bound of the norm of the Fr'echet derivative
of the matrix sign function of A at X given by Roberts [30]:
is the length of the closed contour \Gamma.
Recently, an asymptotic perturbation bound of sign(A) was given by Byers, He
and Mehrmann [13]. They show that to first order in ffi A,
kffiAk;
6 Z. BAI AND J. DEMMEL
where A is assumed to have the form of (1), kffiAk is sufficiently small and
the separation of the matrices A 11 and A 22
[33].\Omega is the Kronecker product. Comparing
the bounds (7) and (9), we note that first the bound (7) is a global bound and
is an asymptotic bound. Second, the assumption (6) for the bound (7) has a simple
geometric interpretation (see Remark 2 above). It is unspecified how to interpret
the assumption on sufficient small kffiAk for the bound (9).
3. Conditioning of matrix sign function computation. In [2], we point out
that it may be much more efficient to compute to half machine precision
only, i.e., to compute S with an absolute error bounded by u 1=2 kSk. To avoid ill
conditioning in the Newton iteration and achieve the half machine precision, we believe
that the matrix A must have condition number less than u \Gamma1=2 . If A is ill conditioned,
say having singular values less than u 1=2 kAk, we need to use a preprocessing step to
deflate small singular values by a unitary similarity transformation, and obtain a
submatrix having condition number less than u \Gamma1=2 , and then compute the matrix
sign function of this submatrix. Such a deflation procedure may be also needed for
the intermediate matrices in the Newton iteration in the worst case.
We now look more closely at the situation of near convergence of the Newton
iteration, and relate the error to the distance to the nearest ill-posed problem [18].
As before, the ill-posed problems are those matrices with pure-imaginary eigenvalues.
Without loss of generality, let us assume A is of the form
A 11 A 12
where Otherwise, for any matrix B, by the Schur
decomposition, we can where A has the above form, and then
R be the solution of the Sylvester equation
A
which must exist and be unique since A 11 and A 22 have no common eigenvalues. Then
it is known that the spectral projector P corresponding to the eigenvalues of A 11 is
' I R
and
. The following lemma relates R and the norm of the projection
P to sign(A) and its condition number.
Lemma 3.1. Let A and R be as above. Let ae
1.
I \Gamma2R
2. ae, and therefore
Proof.
1. Let
' I R
I
. It is easy to verify that if R satisfies (12), then
I \Gamma2R
MATRIX SIGN FUNCTION FOR COMPUTING INVARIANT SUBSPACES 7
2. Using the singular value decomposition (SVD) of R: URV
one can reduce computing the SVD of S to computing the SVD of
I \Gamma2\Sigma
which, by permutations, is equivalent to computing the SVDs of the 2 by 2
matrices
. This is in turn a simple calculation.
We note that for the solution R of the Sylvester equation (12), we have
where the equality is attainable [33]. From Lemma 3.1, we see that the conditioning of
the matrix sign function computation is closely related to the norm of the projection P ,
therefore the norm of R, which in turn is closely related to the quantity
Specifically, when kRk is large,
and
If kA 12 k is moderate, an ill conditioned matrix sign function means large kRk, which
in turn means small Following Stewart [33], it means that it is harder
to separate the invariant subspaces corresponding to the matrices A 11 and A 22 .
The following theorem discusses the conditioning of the eigenvalues of sign(A),
and the distance from sign(A) to the nearest ill-posed problem.
Theorem 3.2. Let A and R be as in Lemma 3.1. Then
1. Let ffi S have the property that S + ffi S has a pure imaginary eigenvalue. Then
may be chosen with no smaller. In the language of
[35], the ffl-pseudospectrum of S excludes the imaginary axis for ffl ! 1=kSk,
and intersects it for ffl - 1=kSk.
2. The condition number of the eigenvalues of S is kPk. In other words, perturbing
S by a small ffi S perturbs the eigenvalues by at most kPk kffiSk+O(kffiSk 2 ).
3. If A is close to S and -(S) ! u \Gamma1=2 , then Newton iteration (2) in floating
point arithmetic will compute S with an absolute error bounded by u 1=2 kSk.
Proof.
1. The problem is to minimize oe min (S \Gamma iiI) over all real i, where oe min is
the smallest singular value of S \Gamma iiI. Using the same unitary similarity
transformation and permutation as in the part 1 of Lemma 3.1, we see that
this is equivalent to minimizing
oe min
over all oe j and real i. This is a straightforward calculation, with the minimum
being obtained for
8 Z. BAI AND J. DEMMEL
2. The condition number of a semi-simple eigenvalue is equal to the secant of the
acute angle between its left and right eigenvectors [24, 17]. Using the above
reduction to 2 by 2 subproblems (this unitary transformation of coordinates
does not changes angles between vectors), this is again a straightforward
calculation.
3. Since the absolute error ffi S in computing 1
essentially by the error in computing S
For the Newton iteration to converge, ffi S cannot be so large that S + ffi S has
pure imaginary eigenvalues; from the result 1, this means kffiSk
Therefore, if u iteration (2)
will compute S with an absolute error bounded by u 1=2 kSk.
It is naturally desired to have an analysis from which we know the conditioning
of the intermediate matrices A k in the Newton iteration. It will help us in addressing
the question of how to detect possible appearance of pure imaginary eigenvalues and
to modify or terminate the iteration early if necessary. Unfortunately, it is difficult to
make a clean analysis far from convergence, because we are unable to relate the error
in each step of the iteration to the conditioning of the problem. We can do a coarse
analysis, however, in the case that the matrix is diagonalizable.
Theorem 3.3. Let A have eigenvalues - j (none pure imaginary or zero), right
eigenvectors x j and left eigenvectors y j , normalized so
Let A k be the matrix obtained at the kth Newton iteration (2). Then for all k,
oe
oe min
Proof. We may express the eigen-decomposition of A as
We
wish to bound j- j;k j from above and below for all k. This is easily done by noting
that
so that all - j;k lie inside a disk defined by
This disk is symmetric about the real axis, so its points of minimum and maximum
absolute value are both real. Solving for these extreme points yields
MATRIX SIGN FUNCTION FOR COMPUTING INVARIANT SUBSPACES 9
This means
oe
Similarly
oe
which proves the bound (15).
As we know, the error introduced at each step of the iteration is mainly caused
by the computation of matrix inverse, which is approximately bounded in norm by
when oe - 1. If uoe \Gamma3 ! oe min error cannot make an intermediate
A k become singular and so cause the iteration to fail. Our analysis shows that if
uoe \Gamma3 ! oe, or oe ? u 1=4 , then the iteration will not fail. This very coarse bound
generalizes result 3 of Theorem 2.
We note that if A is symmetric, by the orthonormal eigendecomposition of
then from Theorem 3, we have
Therefore,
if
It shows that when A is symmetric, the condition number of the intermediate matrices
A k , which affects the numerical stability of the Newton iteration, is essentially
determined by the square of the distance of the eigenvalues to the imaginary axis. 1
When A is nonsymmetric and diagonalizable, from Theorem 3.3, we also see
that the condition number of the intermediate matrices A k is related to the norms
of the spectral projectors
corresponding to the eigenvalues - j
and the quantities of the form
~
by a simple algebraic manipulation,
we have
referee predicted that in the symmetric case, the condition number of A k might be determined
only by the distance, not the square of the distance. We were not able to prove such prediction.
Z. BAI AND J. DEMMEL
From this expression, we see that if there is an eigenvalue - j of A very near to the
pure imaginary axis, i.e, ff j is small, then by the first order Taylor expansion of ~
oe j in
term of ff j , we have
Therefore, to first order in ff j , the condition numbers of the intermediate matrices A k
O
This implies that even if the eigenvalues of A are well-conditioned (i.e., the kP j k are
not too large), if there are also eigenvalues of A closer to the imaginary axis than u 1=2 ,
then the condition number of A k could be large, so the Newton
iteration could fail to converge.
4. Backward Stability of Computed Invariant Subspace. As discussed in
the previous section, because of possible ill conditioning of a matrix with respect to
inversion and rounding errors during the Newton iteration, we generally only expect
to be able to compute the matrix sign function to the square root of the machine
precision, provided that the initial matrix A has condition number smaller than u \Gamma1=2 .
This means that when Newtion iteration converges, the computed matrix sign function
u)kSk:
Under this assumption, b
is an approximate spectral projection corresponding
to -+ (A). Therefore, if
P ), the first ' columns b
in the rank revealing QR decomposition of b
span an approximate invariant
subspace. b
Q has the form
A 11
A 12
A 22
with -( b
being the approximate eigenvalues of A in C+ , and -( b
A 22 ) being the
approximate eigenvalues of A in C \Gamma . Since we expect the computed matrix sign
function to be of half machine precision, it is reasonable to expect computing the
invariant subspace to half precision too. This in turn means that the backward error
in the computed decomposition b
Q is bounded by O(
that the problem is not very ill conditioned. In this section, we will try to justify such
expectation.
To this end, we first need to bound the error in the space spanned by the leading
columns of the transformation matrix Q, i.e., we need to know how much
a right singular subspace of the exact projection matrix perturbed
when P is perturbed by a matrix of norm j. Since P is a projector, the subspace is
spanned by the right singular vectors corresponding to all nonzero singular values of
P (call the set of these singular values S). In practice, of course, this is a question of
rank determination. From the well-known perturbation theory of the singular value
MATRIX SIGN FUNCTION FOR COMPUTING INVARIANT SUBSPACES 11
decomposition [34, page 260], the space spanned by the corresponding singular vectors
is perturbed by at most O(j)=gap S , where gap S is defined by
To compute gap S , we note that there is always a unitary change of basis in which
a projector is of the form
I \Sigma
, where diagonal with
straightforward calculation, we find that the singular values
of the projector are f
where the number
of ones in the set of singular values is equal to maxf2' \Gamma n; 0g. Since
f
Thus, the error ffi Q in Q is bounded by
O(
Hence, the backward error in the computed spectral decomposition is bounded by
is the second order perturbation term of kffiQk. Therefore, if 2' - n, we
have the following first order bound on the backward stability of computed invariant
subspace:
O(
u)kSk
O(
u)kSk
If we use the bound (5) of the matrix sign function S, then from (21), we have
O(
dA
where dA , defined in (4), is the distance to the ill-posed problem. On the other hand,
if we use the bound (13) for the matrix sign function S, then from (21) again, we have
O(
22 ) is the separation of the matrices A 11 and A 22 , if A is assumed
to have the form (11). We note that the error bound (23) is essentially the same as
the error bound given by Byers, He and Mehrmann [13], although we use a different
approach. In [13], it is assumed that in (19), where F 21 is the (2,1)
block of the matrix F . Therefore, O(
u) term in (23) is replaced by O(u).
Z. BAI AND J. DEMMEL
The bounds (22) and (23) reveal two important features of the matrix sign function
based algorithm for computing the invariant subspace. First, they indicate that
the backward error in the computed approximate invariant subspace appears no larger
than the absolute error in the computed matrix sign function, provided that the spectral
decomposition problem is not very ill conditioned (i.e., dA or ffi is not tiny).
Second, if 2' - n, the backward error is a decreasing function of oe l . If oe ' is large, this
means oe 1 and so
are large, and this in turn means the eigenvalues
close the imaginary axis are ill conditioned. It is harder to divide these eigenvalues.
Of course as they become ill conditioned, dA decreases at the same time, which must
counterbalance the increase in oe ' in a certain range.
It is interesting to ask which error bound (22) and (23) is sharper, i.e., which
one of the quantities dA and 22 ) is larger. In [13], an example of a
2 by 2 matrix is given to show that the quantity ffi is larger than the quantity dA .
However, we can also devise simple examples to show that dA can be larger than
More generally, by choosing A 11 to be a large Jordan block with a tiny eigenvalue, and
dA to be close to the square root of ffi . dA is computed using "numerical
brute force" to plot the function dA (- ) on a wide the range of - 2 IR, and search for
the minimal value.
Note that by modifying A to be A \Gamma oeI, where oe is a (sufficiently small) real
number, dA will change, but ffi will not. Thus dA and ffi are not completely comparable
quantities. We believe dA to be a more natural quantity to use than ffi , since ffi does not
always depend on the distance to the nearest ill-posed problem. This is reminiscent
of the difference between the quantities
In practice, we will use the a posteriori bound kE 21 k=kAk anyway, since if we
block upper-triangularize b
Q by setting the (2; 1) block to zero, kE 21 k=kAk is
precisely the backward error we introduce.
Before ending of this section, let us comment on stability of the matrix sign function
based algorithm versus the QR algorithm. The QR algorithm is a numerical
backward stable method for computing the Schur decomposition of a general non-symmetric
matrix A. The computed Schur form b
T and Schur vectors b
Q by the QR
algorithm satisfy
where E is of the order of ukAk. Numerical software for the QR algorithm is available
in EISPACK [32] and LAPACK [1]. Although nonconvergent examples have been
found, they are quite rare in practice [6, 16]. We note that the eigenvalues on the
(block)-diagonal of b
may appear in any order. Therefore, if an application requires
an invariant subspace corresponding to the eigenvalues in a specific region in complex
plane, a second step of reordering eigenvalues on the diagonal of b
T is necessary. A
guaranteed stable implementation of this reordering is described in [7].
The matrix sign function based algorithm can be regarded as an algorithm to
combine these two steps into one. If the matrix sign function can be computed
MATRIX SIGN FUNCTION FOR COMPUTING INVARIANT SUBSPACES 13
within the order of ukSk, then the analysis in this section shows that the matrix
sign function based algorithm could be as stable as the QR algorithm plus reordering.
Unfortunately, if the matrix is ill conditioned with respect to matrix inversion (which
does not affect the QR algorithm), numerical unstable is anticipated in the computed
matrix sign function. Therefore, in general, the matrix sign function is less stable
than the QR algorithm plus reordering.
5. Numerical Experiments. In this section, we will present numerical examples
to verify the above analysis. We will see that the numerical stability of the
Newton iteration (2) and the backward accuracy of computed spectral decomposition
(1) under the influence of the conditioning of the matrix A with respect to inversion,
the condition number -(S) of and the distance \Delta(A) of the eigenvalues
of A to the pure imaginary axis, where We use the easily
computed quantity \Delta(A) as an surrogate of the quantity dA in (4).
Let us recall that the analysis of sections 3 and 4 essentially claims that
(1) If \Delta(A) ! u 1=2 , then the Newton iteration may fail to converge or fail to
compute the matrix sign function within the absolute error u 1=2 kSk, even
when the matrix sign function is well-conditioned. See (18).
(2) If -(S) ? u \Gamma1=2 , then even the distance \Delta(A) is not small, the Newton
iteration may still fail to compute the matrix sign function in the absolute
error of O(u 1=2 kSk). See the part 3 of Theorem 3.2.
(3) In general, the backward error in the computed spectral decomposition will
be smaller than the absolute error in the computed matrix sign function. See
(21).
The following numerical examples will illustrate these claims. Our numerical experiments
were performed on a SUN workstation 10 with machine precision "
u. All the algorithms are implemented in Matlab 4.0a. We use the
simple Newtion iteration (2) to compute the matrix sign function with the stopping
criterion
The maximal number of iterations is set to be 70. At the convergence, we have
S, the computed matrix sign function. We use the QR decomposition
with column pivoting as the rank revealing scheme: 1( b
R\Pi, and finally
compute
A 11
A 12
A 22
where the first
R) columns of b
Q spans the invariant subspaces corresponding
to -( b
A 11 ), which are the approximate eigenvalues of A in C+ . kE 21 k=kAk is the
backward error committed by the algorithm.
All our test matrices are constructed of the form
where U is an orthogonal matrix generated from the QR decomposition a random
matrix with normal distribution having mean 0.0 and variance 1.0. We will choose
different submatrices A 11 , A 22 and A 12 so that the generated matrices A have different
specific features in order to observe our theoretical results in practice.
14 Z. BAI AND J. DEMMEL

Table
Numerical Results for Example 1
The exact matrix sign function of A and the condition number of S
are computed as described in Lemma 3.1. The condition number of A is computed
by Matlab function cond.
In the following tables, iter is the number of iterations of the Newton iteration.
number 10 ff in parenthesis next to an iteration number iter indicates that the
convergence of the Newton iteration was stationary about O(10 ff ) from the iter th
iteration forward, and failed to satisfy the stopping criterion even after the allowed
maximal number of iterations.
We have experimented numerous matrices with different pathologically ill conditioning
in terms of the distance to the pure imaginary axis, the condition numbers
of -(A) and -(S), and the different values of sep(A 11 ; A 22 ) and so on. Two selected
examples presented here are of typical behaviors we observed.
Example 1. In this example, the matrices A are of the form (24) with
and A is a random 2 by 2 matrix with normal distribution
multiplying by a parameter c. The generated matrix A have two complex
conjugate eigenpairs s \Sigma i and \Gammas \Sigma i. As s ! 0, the distance too. The
size of the parameter c will adjust the conditioning of the resulted matrix A and its
matrix sign function.

Table

1 reports the computed results for different values of From the
table, we see that when the matrices are well conditioned and the corresponding the
matrix sign function is also well-conditioned, as stated in the claim (1), the convergence
rate and accuracy of the Newton iteration is clearly determined by the distance
\Delta(A). When the distance becomes smaller, there is a steady increase in the number
of the Newton iteration required to convergence and the loss of the accuracy in the
computed matrix sign function and therefore the desired invariant subspace. From the
table, we also see that when both \Delta(A) and -(S) are moderate, the Newton iteration
fails to compute the matrix sign function in half machine precision. Nevertheless, the
computed invariant subspace seems still have half machine precision, see the claim
MATRIX SIGN FUNCTION FOR COMPUTING INVARIANT SUBSPACES 15

Table
Numerical Results of Example 2
Example 2. In this example, the test matrices A are of the form (24). A 12 are 5 by 5
0)-normally distributed random matrices. The submatrices A 11 and A 22 are first
set by 5 by 5 (1; 0)-normally distributed random upper tridiagonal matrices, and then
the diagonal elements of A 11 and A 22 are replaced by dja ii j and \Gammadja ii j, respectively,
where a ii (1 - i - n) are random numbers with normal distribution (0; 1), d is a
positive parameter. A 12 are 5 by 5 (1; 0)-normally distributed random matrices.
The numerical results are reported in Table 2. For the given parameter d, the
eigenvalues are well-separated away from the pure imaginary axis (\Delta(A) is not small),
however, as stated in the claim (2), we see the influence of the condition numbers -(S)
to the convergence of the Newton iteration and therefore the accuracy of the computed
matrix sign function and the invariant subspace.
6. Refining Estimates of Approximate Invariant Subspaces. When we
use the matrix sign function based algorithm to deflate an invariant subspace of matrix
A, we end up with the form
A 11
A 12
A 22
where the size of kE 21 k=kAk reveals the accuracy and backward stability of computed
invariant subspace spanning by b
of A. If higher accuracy is desired, we may use
iterative refinement techniques to improve the accuracy of computed invariant sub-
space. The methods are due to Stewart [33], Dongarra, Moler and Wilkinson [20], and
Chatelin [15]. Even though these methods all apparently solve different equations,
as shown by Demmel [19], after changing variables, they all solve the same Riccati
equation in the inner loop.
Let us follow Stewart's approach to present the first class of methods. From
(25), we know that b
spans an approximate invariant subspaces and b
spans an
orthogonal complementary subspace. If let the true invariant subspace is represented
as b
its orthogonal complementary subspace as b
Then Y is derived as follows: b
will be an invariant subspace if and only if
the lower left block of
is zero, i.e. if the lower left corner of
I \GammaY H
Y I
A 11
A 12
A 22
\GammaY I
Z. BAI AND J. DEMMEL
is zero. Thus, Y must satisfy the equation
A 12 Y;
which is the well-known algebraic Riccati equation. We may use the following two
iterative methods to solve it:
1. The simple Newton iteration:
A 22 Y
with
2. The modified Newtion iteration:
with
Therefore, we only need to solve a Sylvester equation in the inner loop of the iterative
refinement.
In following numerical example, we only use the simple Newton iteration (26) to
refine the approximate invariant subspace computed by matrix sign function based
algorithm, with the following stopping criterion:
Example 3 We continue the Example 2. Table 3 lists the 22 ), the number of
iterative refinement steps and the backward accuracy of improved invariant subspace.
As shown in the convergence analysis for the iterative solvers (26) and (27) of the
Riccati equation by Stewart [33] and Demmel [18], if we let
A 22 );
then under the assumptions k ! 1=4 and k ! 1=12, the iterations (26) and (27)
converge, respectively. Therefore, sep( b
A 22 ) is a key factor to the convergence of
the iterative refinement schemes. The above examples verify such analysis. From the
analysis of section 3, we recall that sep( b
A 22 ) also affects the backward stability
of the computed invariant subspace by the matrix sign function based algorithm in
the first place (before iterative refinement).
7. Extension to the Generalized Eigenproblem. In this section, we outline
a scheme to extend the matrix sign function based algorithm to solve the generalized
eigenvalue problem of a regular matrix pencil A \Gamma -B. A matrix pencil A \Gamma -B is
regular if A\Gamma-B is square and det(A\Gamma-B) is not identically zero. In [22], Gardiner and
Laub have considered an extension of the Newton iteration for computing the matrix
sign function to a matrix pencil for solving generalized algebraic Riccati equations.
Here we discuss another possible approach, which includes the computation of both
left and right deflating subspaces.
For the given matrix pencil A \Gamma -B, the problem of the spectral decomposition
is to seek a pair of left and right deflating subspaces L and R corresponding to the
eigenvalues of the pencil in a specified region D in complex plane. In other words,
we want to find a pair of unitary matrices QL and QR so that if
MATRIX SIGN FUNCTION FOR COMPUTING INVARIANT SUBSPACES 17

Table
Iterative Refinement Reresults of Example 2
d
28
where the eigenvalues of A 11 \Gamma -B 11 are the eigenvalues of A \Gamma -B in a selected region
D in complex plane. Here, we will only discuss the region D to be the open right half
complex plane. As the same treatment in the standard eigenproblem, by employing
M-obius transformations (ffA divide-and-conquer, D can be
the union of intersections of arbitrary half planes and (complemented) disks, and so
a rather general region.
To this end, by directly applying the Newton iteration to AB \Gamma1 , we have
At convergence, In practice, we do not want to invert B if it is ill
conditioned. Hence, by letting Z then the above iteration becomes
This leads to the following iteration:
converges quadratically to a matrix Z1 . Then
Next, to find the desired deflating subspace, we use
the rank revealing QR decomposition to calculate the range space of the projection
corresponding to the spectral in the open right half plane, which
has the same range space as Thus by computing the rank revealing
QR decomposition of Z1 we obtain the invariant subspace of AB \Gamma1
without inverting B, i.e.,
where -(CR ) are the eigenvalues of the pencil A \Gamma -B in the open right half plane,
-(CL ) are the ones of A \Gamma -B in the open left half plane. Therefore, we have obtained
the left deflating subspace of A \Gamma -B.
Z. BAI AND J. DEMMEL
To compute the right deflating subspace of A\Gamma-B, we can applying the above idea
to A H \Gamma -B H , since transposing swaps right and left spaces. The Newton iteration
implicitly applying to A H B \GammaH turns out to be
~
converges quadratically to a matrix ~
Z1 . Using
the same arguments as above, after computing the rank revealing QR decomposition
of ~
QRRR \Pi R , we have
~
R A H B \GammaH ~
where -(DL ) are the eigenvalues of the pencil A\Gamma-B in the open left half plane, -(DR )
are the ones of A \Gamma -B in the open right half plane. Note that for the desired spectral
decomposition, after transposing, we need to first compute the deflating subspace
corresponding to the eigenvalues in the open left half plane. Let
\Pi, where
~
\Pi is anti-diagonal identity matrix 2 , then we have
From (29) and (30), we immediately have
R D H0 D H
L AQR and Q H
BQR have the partitions
A 21 A 22
we have
R D H0 D H
Note that -(CL ) are the eigenvalues of the pencil A \Gamma -B in the open left half plane,
-(DR ) are the eigenvalues of the pencil A \Gamma -B in the open right half plane. Therefore,
the above homogeneous Sylvester equation has only solution B From (31) or
(32), we have A The computed unitary orthogonal matrices QL and QR give
the desired spectral decomposition (28).
2 The permutation ~
\Pi can be avoided if we use the rank revealing QL decomposition.
MATRIX SIGN FUNCTION FOR COMPUTING INVARIANT SUBSPACES 19
8. Closing Remarks. In this paper, we have presented a number of new results
and approaches to further analyze the numerical behavior of the matrix sign function
and algorithms using it to compute spectral decompositions of nonsymmetric matri-
ces. From this analysis and numerical experiments, we conclude that if the spectral
decomposition problem is not ill conditioned, the algorithm is a practical approach
to solve the nonsymmetric eigenvalue problem. Performance evaluation of the matrix
sign function based algorithm on parallel distributed memory machines, such as the
Intel Delta and CM-5, is reported in [4].
During the course of this work, we have discovered a new approach which essentially
computes the same spectral projection matrix as the matrix sign function
approach does, and also uses basic matrix operations, namely, matrix multiplication
and the QR decomposition. However, it avoids the matrix inverse. From the point of
view of accuracy, this is a more promising approach. The new approach is based on
the work of Bulgakov, Godunov and Malyshev [10, 27, 28]. In [5], we have improved
their results in several important ways, and made it a truly practical and inverse free
highly parallel algorithm for both the standard and generalized spectral decomposition
problems. In brief, the difference between the matrix sign function and inverse
methods is as follows. The matrix sign function method is significantly faster
than it converges, but there are some very difficult problems
where the inverse free algorithm gives a more accurate answer than the matrix
sign function algorithm. The interested reader may see the paper [5] for details.

Acknowledgements

. The first author was supported in part by ARPA grant
DM28E04120 and P-95006 via a subcontract from Argonne National Laboratory and
by an NSF grant ASC-9313958 and in part by an DOE grant DE-FG03-94ER25219 via
subcontracts from University of California at Berkeley. The second author was funded
in part by the ARPA contract DAAH04-95-1-0077 through University of Tennessee
subcontract ORA7453.02, ARPA contract DAAL03-91-C-0047 through University of
Tennessee subcontract ORA4466.02, NSF contracts ASC-9313958 and ASC-9404748,
contracts DE-FG03-94ER25219, DE-FG03-94ER25206, DOE contract No. W-
through subcontract No. 951322401 with Argonne National Labora-
tory, and NSF Infrastructure Grant Nos. CDA-8722788 and CDA-9401156.
The information presented here does not necessarily reflect the position or the
policy of the Government and no official endorsement should be inferred.
The authors would like to acknowledge Ralph Byers, Chunyang He, Nick Higham
and Volker Mehrmann for fruitful discussions on the subject. We would also like to
thank the referees for their valuable comments on the manuscript.



--R


Design of a parallel nonsymmetric eigenroutine toolbox
Design of a parallel nonsymmetric eigenroutine toolbox
The spectral decomposition of nonsymmetric matrices on distributed memory parallel computers.
Inverse free parallel spectral divide and conquer algorithms for nonsymmetric eigenproblems.
Convergence of the shifted QR algorithm on 3 by 3 normal matrices.
Reordering diagonal blocks in real schur form.
A regularity result for the singular values of a transfer matrix and a quadratically convergent algorithm for computing its L1-norm
A bisection method for computing the H1 norm of a transfer matrix and related problems.
Circular dichotomy of the spectrum of a matrix.
Solving the algebraic Riccati equation with the matrix sign function.
A bisection method for measuring the distance of a stable matrix to the unstable matrices.
The matrix sign function method and the computation of invariant subspaces.
On the stability radius of a generalized state-space system
Simultaneous Newton's iteration for the eigenproblem.
How the QR algorithm fails to converge and how to fix it.
The condition number of equivalence transformations that block diagonalize matrix pencils.
On condition numbers and the distance to the nearest ill-posed problem
Three methods for refining estimates of invariant subspaces.
Improving the accuracy of computed eigenvalues and eigenvectors.
Matrix Eigensystem Routines - EISPACK Guide Extension
A generalization of the matrix-sign function solution for algebraic Riccati equations
The matrix sign decomposition and its relation to the polar decomposition.
Perturbation Theory for Linear Operators.
Polar decomposition and matrix sign function condition estimates.
Matrix sign function algorithms for Riccati equa- tions
Guaranteed accuracy in spectral problems of linear algebra
Parallel algorithm for solving some spectral problems of linear algebra.
Condition estimation for the matrix function via the Schur decomposition.
Linear model reduction and solution of the algebraic Riccati equation.
Separation of matrix eigenvalues and structural decomposition of large-scale systems
Matrix Eigensystem Routines - EISPACK Guide
and perturbation bounds for subspaces associated with certain eigenvalue problems.
Matrix Perturbation Theory.
Pseudospectra of matrices.
--TR

--CTR
Daniel Kressner, Block algorithms for reordering standard and generalized Schur forms, ACM Transactions on Mathematical Software (TOMS), v.32 n.4, p.521-532, December 2006
newton's method;matrix sign function;eigenvalue problem;deflating subspaces;invariant subspace