The Matrix Sign Function Method and the Computation of Invariant Subspaces.
A perturbation analysis shows that if a numerically stable procedure is used to compute the matrix sign function, then it is competitive with conventional methods for computing invariant subspaces. Stability analysis of the Newton iteration improves an earlier result of Byers and confirms that ill-conditioned iterates may cause numerical instability. Numerical examples demonstrate the theoretical results.
Introduction
. If A 2 R n\Thetan has no eigenvalue on the imaginary axis, then the
matrix sign function sign(A) may be defined as
Z
(1)
where fl is any simple closed curve in the complex plane enclosing all eigenvalues of A
with positive real part. The sign function is used to compute eigenvalues and invariant
subspaces [2, 4, 6, 13, 14] and to solve Riccati and Sylvester equations [9, 15, 16, 28]. The
matrix sign function is attractive for machine computation, because it can be efficiently
evaluated by relatively simple numerical methods. Some of these are surveyed in [28].
It is particularly attractive for large dense problems to be solved on computers with
advanced architectures [2, 11, 16, 33].
Beavers and Denman use the following equivalent definition [6, 13]. Let
be the Jordan canonical decomposition of a matrix A having no eigenvalues
on the imaginary axis. Let the diagonal part of J be given by the matrix
then
(A) be the invariant subspace of A corresponding to eigenvalues
with positive real part, let be the invariant subspace of A corresponding
to eigenvalues with negative real part, let be the skew projection onto
parallel to be the skew projection onto parallel to
. Using the same contour fl as in (1), the projection P + has the resolvent integral
representation [23, Page 67] [2]
Z
(2)
To appear in SIAM Journal on Matrix Analysis and Its Applications.
y University of Kansas, Dept. of Mathematics, Lawrence, Kansas 66045, USA. Partial support was
received from National Science Foundation grants INT-8922444 and CCR-9404425 and University of
Kansas GRF Allocation 3514-20-0038.
z TU-Chemnitz-Zwickau, Fak. f. Mathematik, D-09107 Chemnitz, FRG. Partial support was received
from Deutsche Forschungsgemeinschaft, Projekt La 767/3-2.
It follows from (1) and (2) that
The matrix sign function was introduced using definition (1) by Roberts in a 1971
technical report [34] which was not published until 1980 [35]. Kato [23, Page 67] reports
that the resolvent integral (2) goes back to 1946 [12] and 1949 [21, 22].
There is some concern about the numerical stability of numerical methods based
upon the matrix sign function [2, 8, 19]. In this paper, we demonstrate that evaluating
the matrix sign function is a more ill-conditioned computational problem than the
problem of finding bases of the invariant subspaces V
Section 3. Nevertheless, we also give perturbation and error analyses, which show
that (at least for Newton's method for the computation of the matrix sign function
[8, 9]) in most circumstances the accuracy is competitive with conventional methods
for computing invariant subspaces. Our analysis improves some of the perturbation
bounds in [2, 8, 18, 24].
In Section 2 we establish some notation and clarify the relationship between the
matrix sign function and the Schur decomposition. The next two sections give a perturbation
analysis of the matrix sign function and its invariant subspaces. Section 5 gives
a posteriori bounds on the forward and backward error associated with a corrupted
value of sign(S). Section 6 is a stability analysis of the Newton iteration.
Throughout the paper, k \Delta k denotes the spectral norm, k \Delta k 1 the 1-norm (or column
sum norm), and k \Delta k F the Frobenius norm k \Delta k
. The set of eigenvalues
of a matrix A is denoted by -(A). The open left half plane is denoted by C \Gamma and the
open right half plane is denoted by C Borrowing some terminology from engineering,
we refer to the invariant subspace of a matrix A 2 R n\Thetan corresponding
to eigenvalues in C \Gamma as the stable invariant subspace and the subspace
corresponding to eigenvalues in C + as the unstable invariant subspace. We use
for the the skew projection onto V + parallel to V \Gamma and for the
skew projection onto
2. Relationship with the Schur Decomposition. Suppose that A has the
Schur form
k A 11 A 12
where Y is a solution of the Sylvester equation
then
and
The solution of (4) has the integral representation
Y =-i
Z
where fl is a closed contour containing all eigenvalues of A with positive real part
[29, 36]).
The stable invariant subspace of A is the range (or column space) of
is a QR factorization with column pivoting [1, 17], where Q and R are partitioned in
the obvious way, then the columns of Q 1 form an orthonormal basis of this subspace.
Here Q is orthogonal, \Pi is a permutation matrix, R is upper triangular, and R 1 is
nonsingular.
It is not difficult to use the singular value decomposition of Y to show that
s
It follows from (4) that
where sep is defined as in [17] by
3. The Effect of Backward Errors. In this section we discuss the sensitivity
of the matrix sign function subject to perturbations. For a perturbation matrix E, we
give first order estimates for E) in terms of submatrices and powers of kEk.
Based on Fr'echet derivatives, Kenney and Laub [24] presented a first order perturbation
theory for the matrix sign function via the solution of a Sylvester equation.
Mathias [30] derived an expression for the Fr'echet derivative using the Schur decom-
position. Kato's encyclopedic monograph [23] includes an extensive study of series
representations and of perturbation bounds for eigenprojections. In this section we
derive an expression for the Fr'echet derivative using integral formulas.
Let
oe min
where oe min is the smallest singular value of A \Gamma -iI. The quantity dA is the
distance from A to the nearest complex matrix with an eigenvalue on the imaginary
axis. Practical numerical techniques for calculating dA appear in [7, 10]. If
then E is too small to perturb an eigenvalue of A on or across the imaginary axis. It
follows that for kEk ! dA , sign(A+E), and the stable and unstable invariant subspaces
of are smooth functions of E.
Consider the relatively simple case in which A is block diagonal.
Lemma 3.1. Suppose A is block diagonal,
where . Partition the perturbation E 2 R n\Thetan conformally
with A as
where F 12 and F 21 satisfy the Sylvester equations
A 22 F
Proof. The hypothesis that implies that the eigenvalues of A 11
have negative real part and the eigenvalues of A 22 positive real part. In
the definition (1) choose the contour fl to enclose neither
In particular, for all complex numbers z lying on the contour
E) are nonsingular and
E) =-i
Z
Z
e
Z
Partitioning F conformally with E and A, then we have
F 11 =2-i
Z
Z
Z
F 22 =2-i
Z
As in (6), F 12 and F 21 are the solutions to the Sylvester equations (11) and (12) [29, 36].
The contour fl encloses no eigenvalues of A 11 , so
inside fl and F
We first prove that F 22 = 0 in the case that A 22 is diagonalizable, say A
where
Z
Each component of the above integral is of the form R
constant c. If then this is the integral of a residue free holomorphic function and
hence it vanishes. If j 6= k, then
Z
c
Z
c
The general case follows by taking limits of the diagonalizable case and using the dominated
convergence theorem.
The following theorem gives the general case.
Theorem 3.2. Let the Schur form of A be given by (3) and partition E conformally
as
and Y satisfies (4), then
e
~
~
satisfies the Sylvester equation
A 22
~
and ~
~
~
Proof. If
A 11 A 12
and
It follows from Lemma 3.1 that
I
~
on the left side and SQ H on the
right side of the above equation, we have
Y ~
It is easy to verify that
Y ~
\GammaI Y
I
~
\GammaI Y
I
The theorem follows from
\GammaI Y
If dA is small relative to kAk or kY k is large relative to kAk, then the hypothesis
that kSES may be restrictive. However, a small value of dA indicates that A is
very near a discontinuity in sign(A). A large value of kY k indicates that
is small and the stable invariant subspace is ill-conditioned [37, 39].
Of course Theorem 3.2 also gives first order perturbations for the projections
Corollary 3.3. Let the Schur form of A be given as in (3) and let E be as in
(13). Under the hypothesis of Theorem 3.2, the projections P
E)
are as in the statement of Theorem 3.2.
Taking norms in Theorem 3.2 gives the first order perturbation bounds of the next
corollary.
Corollary 3.4. Let the Schur form of A be given as in (3), E as in (13) and
then the first order perturbation of the matrix sign function
stated in Theorem 3.2 is bounded by
On first examination, Corollary 3.4 is discouraging. It suggests that calculating the
matrix sign function may be more ill-conditioned than finding bases of the stable and
unstable invariant subspace. If the matrix A whose Schur decomposition appears in
(3) is perturbed to A+ E, then the stable invariant subspace, range(Q 1 ), is perturbed
to 3.4 and the following
example show that sign(A+E) may indeed differ from sign(A) by a factor of ffi \Gamma3 which
may be much larger than kEk=ffi.
Example 1. Let
The matrix A is already in Schur form, so
have
E) =p
The difference is
Perturbing A to A+E does indeed perturb the matrix sign function by a factor of ffi \Gamma3 .
Of course there is no rounding error in Example 1, so the stable invariant subspace
of A+E is also the stable invariant subspace of sign(A+E) and, in particular, evaluating
exactly has done no more damage than perturbing A. The stable invariant
subspace of A is
"0
); the stable invariant subspace of A
E) is
For a general small perturbation matrix E, the angle between
is of order no larger than O(1=j) [17, 37, 39]. The matrix sign function (and the
projections may be significantly more ill-conditioned than the stable and
unstable invariant subspaces. Nevertheless, we argue in this paper that despite the
possible poor conditioning of the matrix sign function, the invariant subspaces are
usually preserved about as accurately as their native conditioning permits.
However, if the perturbation E is large enough to perturb an eigenvalue across or
on the imaginary axis, then the stable and unstable invariant subspaces may become
confused and cannot be extracted from E). This may occur even when the
invariant subspaces are well-conditioned, since the sign function is not defined in this
case. In geometric terms, in this situation A is within distance kEk of a matrix with
an eigenvalue with zero real part. This is a fundamental difficulty of any method that
identifies the two invariant subspaces by the signs of the real parts of the corresponding
eigenvalues of A.
4. Perturbation Theory for Invariant Subspaces of the Matrix Sign Func-
tion. In this section we discuss the accuracy of the computation of the stable invariant
subspace of A via the matrix sign function.
An easy first observation is that if the computed value of sign(A) is the exact value
of sign(A+E) for some perturbation matrix E, then the exact stable invariant subspace
of E) is also an invariant subspace of A+ E. Let A have Schur form (3) and
let E be a perturbation matrix partitioned conformally as in (13). Let Q 1 be the first
k columns of Q and Q 2 be the remaining columns. If
then A has stable invariant subspace has an invariant
subspace
where ffl 39]. The singular values of W are the tangents of the
canonical angles between In particular, the
canonical angles are at most of order O(1=
Unfortunately, in general, we cannot apply backward error analysis, i.e. we cannot
guarantee that the computed value of sign(A) is exactly the value of E) for
some perturbation E. Consider instead the effect of forward errors, let
where F represents the forward error in evaluating the matrix sign function. Let A have
Schur form (3). Partition Q H sign(A)Q and Q H FQ as
and
where Q is the unitary factor from the Schur decomposition of A (3) and Y is a solution
of (4).
Assume that
and let OE
Perturbing sign(A) to changes the invariant subspace
and by (8) and
s
A
- OE 21
Since obeys the bounds
- 4OE 21
Comparing (19) with (20) we see that the error bound (20) is no greater than twice
the error bound (19). Loosely speaking, a small relative error in sign(A) of size ffl might
perturb the stable invariant subspace by not much more than twice as much as a relative
error of size ffl in A can.
Therefore, the stable and unstable invariant subspaces of sign(A) may be less ill-conditioned
and are never significantly more ill-conditioned than the corresponding
invariant subspaces of A. There is no fundamental numerical instability in evaluating
the matrix sign function as a means of extracting invariant subspaces. However, numerical
methods used to evaluate the matrix sign function may or may not be numerically
unstable.
Example 1 continued. To illustrate the results, we give a comparison of our
perturbation bounds and the bounds given in [3] for both the matrix sign function and
the invariant subspaces in the case of Example 1. The distance to the nearest ill-posed
problem, i.e., is the smallest singular
value of leads to an overestimation of the error in [3]. Since dA - j \Gamma2 , the
bounds given in [3] are, respectively, O(j \Gamma4 ) for the matrix sign function and O(j \Gamma2 )
for the invariant subspaces.
5. A Posteriori Backward and Forward Error Bounds. A priori backward
and forward error bounds for evaluation of the matrix sign function remain elusive.
However, it is not difficult to derive a posteriori error bounds for both backward and
forward error.
We will need the following lemma to estimate the distance between a matrix S and
sign(S).
Lemma 5.1. If S 2 R n\Thetan has no eigenvalue with zero real part and k sign(S)S
Proof. Let S. The matrices F , S, and sign(S) commute, so
This implies that
Taking norms and using kFS
and the lemma follows.
It is clear from the proof of the Lemma 5.1 that
asymptotically correct as k sign(S) \Gamma Sk tends to zero. The bound in the lemma tends
to overestimate smaller values of k sign(S) \Gamma Sk by a factor of two.
Suppose that a numerical procedure for evaluating sign(A) applied to a matrix
A 2 R n\Thetan produces an approximation S 2 R n\Thetan . Consider the problem of finding
small norm solutions E 2 R n\Thetan and F 2 R n\Thetan to Of course, this
does not uniquely determine E and F . Common algorithms for evaluating sign(A) like
Newton's method for the square root of I guarantee that S is very nearly a square root
of I [19], i.e., S is a close approximation of sign(S). In the following theorem, we have
arbitrarily taken
Theorem 5.2. If k sign(S)S
perturbation matrices E and F satisfying
and
Proof. The matrices S +F and A+E commute, so an underdetermined, consistent
system of equations for E in terms of S, A, and
E(S
Let
\GammaI Y
I
be a Schur decomposition of sign(S) whose unitary factor is U and whose triangular
factor is on the right-hand-side of (24). Partition U H EU and U H AU conformally with
the right-hand-side of (24) as
and
A 11 A 12
A 21 A 22
Multiplying (23) on the left by U H and on the right by U and partitioning gives
Y A 21 \GammaA
One of the infinitely many solutions for E is given by
For this choice of E, we have
from which the theorem follows.
Lemma 5.1 and Theorem 5.2 agree well with intuition. In order to assure small
forward error, S must be a good approximate square root of I and, in addition, to assure
small backward error, sign(S) must nearly commute with the original data matrix A.
Newton's method for a root of I tends to do a good job of both [19]. (Note
that in general, Newton's method makes a poor algorithm to find a square root of a
matrix. The square root of I is a special case. See [19] for details.) In particular, the
hypothesis that k sign(S)S usually satisfied when the matrix sign function
is computed by the Newton algorithm.
When S - sign(S), the quantity k
S+S \Gamma1j
S+S \Gamma1j
k makes a good estimate
of the right-hand-side of (21). The bound (22) is easily computed or estimated from
the known values of A and S. However, these expressions are prone to subtractive
cancellation of significant digits.
The quantity kE 21 k is related by (18) to perturbations in the stable invariant sub-
space. The bounds (21) and (22) are a fortiori bounds on kE 21 k, but, as the (1; 2) block
of (25) suggests, they tend to be pessimistic overestimates of kE 21 k if kSk AE 1.
6. The Newton Iteration for the Computation of the Matrix Sign Func-
tion. There are several numerical methods for computing the matrix sign function
[2, 25]. Among the simplest and most commonly used is the Newton-Raphson method
for a root of starting with initial guess It is easily implemented
using matrix inversion subroutines from widely available, high quality linear algebra
packages like LAPACK [1, 2]. It has been extensively studied and many variations have
been suggested [2, 4, 5, 9, 18, 27, 25, 26, 28].
Algorithm 1. Newton Iteration (without scaling)
If A has no eigenvalues on the imaginary axis, then Algorithm 1 converges globally
and locally quadratically in a neighborhood of sign(A) [28]. Although the iteration
ultimately converges rapidly, initially convergence may be slow. However, the initial
convergence rate (and numerical stability) may be improved by scaling [2, 5, 9, 18, 27,
25, 26, 28]. A common choice is to scale X k 1=j det(X k )j (1=n) [9].
Theorem 3.2 shows that the first order perturbation of sign(A) may be as large as
is the relative uncertainty in A. (If there is no other uncertainty,
then ffl is at least as large as the round-off unit of the finite precision arithmetic.) Thus,
it is reasonable to stop the Newton iteration when
The ad hoc constant C is chosen in order to avoid extreme situations, e.g.,
This choice of C works well in our numerical experiments up to
shows furthermore that it is often advantageous to take an extra step of the iteration
after the stopping criterion is satisfied.
Example 2. This example demonstrates our stopping criterion. Algorithm 1 was
implemented in MATLAB 4.1 on an HP 715/33 workstation with floating point relative
accuracy We constructed a
is a random unitary matrix and R an upper triangular matrix with diagonal elements
parameter ff in the (k; k
position and zero everywhere else. We chose ff such that the norm k sign(A)k 1 varies
from small to large.
The typical behavior of the error is that it goes down and then
becomes stationary. This behavior is shown in the Figure 1 for the cases
and
Stopping criterion (26) is satisfied with 1000n at the 8-th step for
at the 7-th step for Taking one extra step would stop at the 9-th
step for and at the 8-th step for
In exact arithmetic, the stable and unstable invariant subspaces of the iterates
are the same as those of A. However, in finite precision arithmetic, rounding errors
perturb these subspaces. The numerical stability of the Newton iteration for computing
the stable invariant subspace has been analyzed in [8], we give an improved error bound
here.
Let X and X be, respectively, the computed k-th and 1)-st iterate of the
Newton iteration starting from
A 11 A 12
Fig. 1. 2.
the number of iterations
log10(||X_k-S||_1)
alpha=0
Suppose that X and X + have the form
A successful rounding error analysis must establish the relationship between E +and
. To do so we assume that some stable algorithm is applied to compute the inverse
in the Newton iteration. More precisely we assume that X
where
for some constant c. Note that this is a nontrivial assumption. Ordinarily, if Gaussian
elimination with partial pivoting is used to compute the inverse, the above error bound
can be shown to hold only for each column separately [8, 38]. The better inversion
algorithms applied to "typical" matrices satisfy this assumption [38, p. 151], but it is
difficult to determine if this is always the case [31, pp. 22-26], [20, p. 150].
Write EX and EZ as
The following theorem bounds indirectly the perturbation in the stable invariant
subspace.
Theorem 6.1. Let X, and EZ be as in (27), (28), (31), and (32). If2
where c is as in (29) and (30), then
Proof. We start with (28). In fact the relationship between E 21 and
from applying the explicit formula for the inverse of (X
~
c
~
c
c
c
Here,
~
Then
~
c
~
c
Using the Neumann lemma that if kBk ! 1, then
have
The following inequalities are established similarly.
22 k
22 k(kE
22 k
Inserting these inequalities in (33) we obtain
The bound in Theorem 6.1 is stronger than the bound of Byers in [8]. It follows
from (19) and Theorem 6.1, that if
then rounding errors in a step of Newton corrupt the stable invariant subspace by no
more than one might expect from the perturbation E 21 in (27). The term sep(X
22
is dominated by
sep
So to guarantee that rounding errors in the Newton iteration do little damage, the
factors in the bound of Theorem 6.1, kX \Gamma1
22 k and (kXk should be
small enough so that
Very roughly speaking, to have numerical stability throughout the algorithm, neither
22 k nor (kXk should be much larger than 1=
The following example from [4] demonstrates numerical instability that can be
traced to the violation of inequality (34).
Example 3. Let
A 11 =6 6 6 6 4
be a real matrix, and let A 22 = \GammaA T. Form
A 11 A 12
and
where the orthogonal matrix Q is chosen to be the unitary factor of the QR factorization
of a matrix with entries chosen randomly uniformly distributed in the interval [0; 1].
The parameter ff is taken as so that there are two eigenvalues of A
close to the imaginary axis from the left and right side. The entries of A 12 are chosen
randomly uniformly distributed in the interval [0; 1], too. The entries of E 21 are chosen
randomly uniformly distributed in the interval [0; eps], where \Gamma16 is the
machine precision.

Table
Evolution of kE
in Example 3.
A R
3 1.2093e-07 2.5501e-08 1.6948e+00
5 7.3034e-08 5.4025e-09 2.0000
6 7.3164e-08 2.7012e-09 2.0000
7 7.2020e-08 1.3506e-09 2.0000
8 7.1731e-08 6.7532e-10 2.0000
9 7.1866e-08 3.3766e-10 2.0000
13 7.1934e-08 2.1151e-11 2.0000
14 7.1938e-08 1.0646e-11 2.0000
19 7.1937e-08 1.7291e-12 2.0000
In this example,
shows the evolution of kE during the Newton iteration starting with
respectively, where E 21 is as in (27). The norm is taken to be
the 1-norm. Because kA \Gamma1k 1
is violated in the first step of the Newton iteration for the starting matrix A, which is
shown in the first column of the table. Newton's method never recovers from this.
It is remarkable, however, that Newton's method applied to R directly seems to
recover from the loss in accuracy in the first step. The second column shows that
although \Gamma7 at the first step, it is reduced by the factor
1=2 every step until it reaches 1:7 \Theta 10 \Gamma12 which is approximately
Observe that in this case the perturbation E 00
in EZ as in (28) is zero and
dominated by 1(kE
22
It is surprising to see that from the second
step on kX \Gamma1
22 k 1 is as small as eps, since A \Gamma1
11 and A \Gamma1
22 do not explicitly appear
in the term X \Gamma1
22 after the first step.
Our analysis suggests that the Newton iteration may be unstable when X k is ill-
conditioned. To overcome this difficulty the Newton iteration may be carried out with
a shift along the imaginary line. In this case we have to use complex arithmetic.
Algorithm 2. Newton Iteration With Shift
END
The real parameter fi is chosen such that oe min fiiI) is not small. For Example
2, when fi is taken to be 0:8, we have
Then by our analysis the computed invariant subspace is guaranteed to be accurate.
7. Conclusions. We have given a first order perturbation theory for the matrix
sign function and an error analysis for Newton's method to compute it. This analysis
suggests that computing the stable (or unstable) invariant subspace of a matrix with
the Newton iteration in most circumstances yields results as good as those obtained
from the Schur form.
8.

Acknowledgments

. The authors would like to express their thanks to N.
Higham for valuable comments on an earlier draft of the paper and Z. Bai and P.
Benner for helpful discussions.



--R


Design of a parallel nonsymmetric eigenroutine toolbox
Design of a parallel nonsymmetric eigenroutine toolbox
Inverse free parallel spectral divide and conquer algorithms for nonsymmetric eigenproblems.
Accelerated convergence of the matrix sign function method of solving Lyapunov
A computational method for eigenvalues and eigenvectors of a matrix with real eigenvalues.
A regularity result and a quadratically convergent algorithm for computing its L1 norm.
Numerical stability and instability in matrix sign function based algorithms.
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.
A systolic algorithm for Riccati and Lyapunov equations.
Perturbations des transformations autoadjointes dans l'espace de Hilbert.
The matrix sign function and computations in systems.
Spectral decomposition of a matrix using the generalized sign matrix.
A generalization of the matrix-sign-function solution for algebraic Riccati equations
Parallel algorithms for algebraic Riccati equations.
Matrix Computations.
Computing the polar decomposition - with applications
Newton's method for the matrix square root.
A survey of error analysis.
On the convergence of the perturbation method
On the convergence of the perturbation method
Perturbation Theory for Linear Operators.
Polar decompositions and matrix sign function condition estimates.
Rational iterative methods for the matrix sign function.

On scaling Newton's method for polar decompositions and the matrix sign function.
The matrix sign function.
The Theory of Matrices.
Condition estimation for the matrix sign function via the Schur decomposition.
Software for Roundoff Analysis of Matrix Algorithms.
Schur complement and statistics.
A parallel algorithm for the matrix sign function.
Linear model reduction and solution of algebraic Riccati equation by use of the sign function.
Linear model reduction and solution of the algebraic Riccati equation by use of the sign function.
and perturbation bounds for subspaces associated with certain eigenvalue problems.
Introduction to Matrix Computations.
Matrix Perturbation Theory.
--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
Peter Benner , Maribel Castillo , Enrique S. Quintana-Ort , Vicente Hernndez, Parallel Partial Stabilizing Algorithms for Large Linear Control Systems, The Journal of Supercomputing, v.15 n.2, p.193-206, Feb.1.2000
matrix sign function;invariant subspaces;perturbation theory