Staircase Failures Explained by Orthogonal Versal Forms.
Treating matrices as points in n2-dimensional space, we apply geometry to study and explain algorithms for the numerical determination of the Jordan structure of a matrix. Traditional notions such as sensitivity of subspaces are replaced with angles between tangent spaces of manifolds in n2-dimensional space. We show that the subspace sensitivity is associated with a small angle between complementary subspaces of a tangent space on a manifold in n2-dimensional space. We further show that staircase algorithm failure is related to a small angle between what we call staircase invariant space and this tangent space. The matrix notions in n2-dimensional space are generalized to pencils in 2mn-dimensional space. We apply our theory to special examples studied by Boley, Demmel, and Kgstrm.
Introduction
. The problem of accurately computing Jordan and Kronecker canonical structures
of matrices and pencils has captured the attention of many specialists in numerical linear algebra.
Standard algorithms for this process are denoted "staircase algorithms" because of the shape of the
resulting matrices [22, Page 370], but understanding of how and why they fail is incomplete. In this
paper, we study the geometry of matrices in n 2 dimensional space and pencils in 2mn dimensional
space to explain these failures. This follows a geometrical program to complement and perhaps replace
traditional numerical concepts associated with matrix subspaces that are usually viewed in n dimensional
space.
This paper targets expert readers who are already familiar with the staircase algorithm. We refer
readers to [22, Page 370] and [10] for excellent background material and we list other literature in Section
1.1 for the reader wishing a comprehensive understanding of the algorithm. On the mathematical side,
it is also helpful if the reader has some knowledge of Arnold's theory of versal forms, though a dedicated
reader should be able to read this paper without such knowledge, perhaps skipping Section 3.2.
The most important contributions of this paper may be summarized:
ffl A geometrical explanation of staircase algorithm failures
Department of Mathematics Room 2-380, Massachusetts Institute of Technology, Cambridge, MA
02139-4307, edelman@math.mit.edu, http://www-math.mit.edu/~edelman, supported by NSF grants 9501278-DMS and
9404326-CCR.
y Department of Mathematics Room 2-333, Massachusetts Institute of Technology, Cambridge, MA 02139-4307,
http://www-math.mit.edu/~yanyuan, supported by NSF grants 9501278-DMS.
ffl Identification of three significant subspaces that decompose matrix or pencil space: T b ,R, S. The
most important of these spaces is S, which we choose to call the "staircase invariant space".
ffl The idea that the staircase algorithm computes an Arnold normal form that is numerically more
appropriate than Arnold's ``matrices depending on parameters''.
ffl A first order perturbation theory for the staircase algorithm
ffl Illustration of the theory using an example by Boley [3]
The paper is organized as follows: In Section 1.1 we briefly review the literature on staircase algorithms.
In Section 1.2 we introduce concepts that we call pure, greedy and directed staircase to emphasize
subtle distinctions on how the algorithm might be used. Section 1.3 contains some important messages
that result from the theory to follow.
Section 2 presents two similar looking matrices with very different staircase behavior. Section 3
studies the relevant n 2 dimensional geometry of matrix space while Section 4 applies this theory to the
staircase algorithm. The main result may be found in Theorem 6.
Sections 5, 6 and 7 mimic Sections 2, 3 and 4 for matrix pencils. Section 8 applies the theory towards
special cases introduced by Boley [3] and Demmel and B. Kagstrom [12].
1.1. Jordan/Kronecker Algorithm History. The first staircase algorithm was given by Kubla-
novskaya for Jordan structure in 1966 [31], where a normalized QR factorization is used for rank determination
and nullspace separation. Ruhe [34] first introduced the use of the SVD into the algorithm in
1970. The SVD idea is further developed by Golub and Wilkinson [23, Section 10]. Kagstrom and Ruhe
[27, 28] wrote the first library quality software for the complete JNF reduction, with the capability of
returning after different steps in the reduction. Recently, Chatitin-Chatelin and Frayss'e [6] developed a
non-staircase "qualitative" approach.
The staircase algorithm for the Kronecker structure of pencils is given by Van Dooren [13, 14, 15] and
Kagstrom and Ruhe [29]. Kublanovskaya [32] fully analyzed the AB algorithm, however, earlier work on
the AB algorithm goes back to the 1970s. Kagstrom [25, 26] gave a RGDSVD/RGQZD algorithm and
this provided a base for later work on software. Error bounds for this algorithm are given by Demmel
and Kagstrom [8, 9]. Beelen and Van Dooren [2] gave an improved algorithm which requires O(m 2 n)
operations for m \Theta n pencils. Boley [3] studied the sensitivity of the algebraic structure. Error bounds
are given by Demmel and Kagstrom [10, 11].
Staircase algorithms are used both theoretically and practically. Elmroth and Kagstrom [19] use
the staircase algorithm to test the set of 2-by-3 pencils hence to analyze the algorithm, Demmel and
Edelman [7] use the algorithm to calculate the dimension of matrices and pencils with a given form.
Dooren [14, 20, 30, 5], Emami-Naeini [20], Kautsky and Nichols [30], Boley [5], Wicks and DeCarlo
[35] consider systems and control applications. Software for control theory is provided by Demmel and
Kagstrom [12].
A number of papers use geometry to understand Jordan and Kronecker structure problems. Fair-
grieve [21] regularizes by taking the most degenerate matrix in a neighborhood, Edelman, Elmroth and
Kagstrom [17, 18] study versality and stratifications, and Boley [4] concentrates on stratifications.
1.2. The Staircase Algorithms. Staircase algorithms for the Jordan and Kronecker form work by
making sequences of rank decisions in combination with eigenvalue computations. We wish to emphasize
a few variations on how the algorithm might be used by coining the terms pure staircase, greedy
staircase, and directed staircase. Pseudocode for the Jordan versions appear near the end of this
subsection. In combination with these three choices, one can choose an option of zeroing or not. These
choices are explained below.
The three variations for purposes of discussion are considered in exact arithmetic. The pure version
is the pure mathematician's algorithm: it gives precisely the Jordan structure of a given matrix. The
greedy version (also useful for a pure mathematician!) attempts to find the most "interesting" Jordan
structure near the given matrix. The directed staircase attempts to find a nearby matrix with a
preconceived Jordan structure. Roughly speaking, the difference between pure, greedy, and directed is
whether the Jordan structure is determined by the matrix, a user controlled neighborhood of the matrix,
or directly by the user respectively.
In the pure staircase algorithm, rank decisions are made using the singular value decomposition. An
explicit distinction is made between zero singular values and nonzero singular values. This determines
the exact Jordan form of the input matrix.
The greedy staircase algorithm attempts to find the most interesting Jordan structure nearby the
given matrix. Here the word "interesting" (or degenerate) is used in the sense of precious gems, the
rarer, the more interesting. Algorithmically, as many singular values as possible are thresholded to zero
with a user defined threshold. The more singular values that are set to 0, the rarer in the sense of
codimension (see [7, 17, 18]).
The directed staircase algorithm allows the user to decide in advance what Jordan structure is
desired. The Jordan structure dictates which singular values are set to 0. Directed staircase is used in a
few special circumstances. For example, it is used when separating the zero Jordan structure from the
right singular structure (used in GUPTRI [10, 11]). Moreover, Elmroth and Kagstrom imposed structures
by the staircase algorithm in their investigation of the set of 2 \Theta 3 pencils [19]. Recently, Lippert and
Edelman [33] use directed staircase to compute an initial guess for a Newton minimization approach to
computing the nearest matrix with a given form in the Frobenius norm.
In the greedy and directed modes if we explicitly zero the singular values, we end up computing a
new matrix in staircase form that has the same Jordan structure as a matrix near the original one. If we
do not explicitly zero the singular values, we end up computing a matrix that is orthogonally similar
to the original one (in the absence of roundoff errors), that is nearly in staircase form. For example,
in GUPTRI [11], the choice of whether to zero the singular values is made by the user with an input
parameter named zero which may be true or false.
To summarize the many choices associated with a staircase algorithm, there are really five distinct
algorithms worth considering: the pure algorithm stands on its own, otherwise the two choices of combinatorial
structure (greedy and directed) may be paired with the choice to zero or not. Thereby we have
the five algorithms:
1. pure staircase
2. greedy staircase with zeroing
3. greedy staircase without zeroing
4. directed staircase with zeroing
5. directed staircase without zeroing
Notice that in the pure staircase, we do not specify zeroing or not, since both will give the same
result vacuously.
Of course algorithms run in finite precision. One further detail is that there is some freedom in
the singular value calculations which lead to an ambiguity in the staircase form: in the case of unequal
singular values, an order must be specified, and when singular values are equal, there is a choice of basis
to be made. We will not specify any order for the SVD, except that all singular values considered to be
zero appear first.
In the ith loop iteration, we use w i to denote the number of singular values that are considered to
be 0. For the directed algorithm, w i are input, otherwise, w i are computed. In pseudocode, we have the
following staircase algorithms for computing the Jordan form corresponding to eigenvalue .
INPUT:
specify pure, greedy, or direct mode
specify zeroing or not zeroing
OUTPUT:
matrix A that may or may not be in staircase form
while A tmp not full rank
Use the SVD to compute an n tmp by n tmp unitary matrix V whose leading w i columns
span the nullspace or an approximation
Choice I: Pure: Use the SVD algorithm to compute w i and the exact nullspace
Choice II: Greedy: Use the SVD algorithm and threshold the small singular values with
a user specified tolerance, thereby defining w i . The corresponding singular vectors
become the first w i vectors of V .
Choice III: Directed: Use the SVD algorithm, the w i are defined from the input Jordan
structure. The w i singular vectors are the first w i columns of V .
Let A tmp be the lower right n corner of A
endwhile
If zeroing, return A in the form I + a block strictly upper triangular matrix.
While the staircase algorithm often works very well, it has been known to fail. We can say that
the greedy algorithm fails if it does not detect a matrix with the least generic form [7] possible within
a given tolerance. We say that the directed algorithm fails if the staircase form it produces is very far
(orders of magnitude, in terms of the usual Frobenious norm of matrix space) from the staircase form
of the nearest matrix with the intended structure. In this paper, we mainly concentrate on the greedy
staircase algorithm and its failure, but the theory is applicable to both approaches. We emphasize that
we are intentionally vague about how "far" is "far" as this may be application dependent, but we will
consider several orders of magnitude to constitute the notion of "far".
1.3. Geometry of Staircase and Arnold forms. Our geometrical approach is inspired by
Arnold's theory of versality [1]. For readers already familiar with Arnold's theory, we point out that we
have a new normal form that enjoys the same properties as Arnold's original form, but is more useful
numerically. For numerical analysts, we point out that these ideas are important for understanding the
staircase algorithm. Perhaps it is safe to say that numerical analysts have had an "Arnold Normal Form"
for years, but we did not recognize as such - the computer was doing it for us automatically.
The power of the normal form that we introduce in Section 3 is that it provides a first order rounding
theory of the staircase algorithm. We will show that instead of decomposing the perturbation space into
the normal space and a tangent space at a matrix A, the algorithm chooses a so called staircase invariant
space to take the place of the normal space. When some directions in the staircase invariant space are
very close to the tangent space, the algorithm can fail.
From the theory, we decompose the matrix space into three subspaces that we call T b , R and S, the
precise definitions of the three spaces are given in Definitions 1 and 3. Here, T b and R are two subspaces
of the tangent space, and S is a certain complimentary space of the tangent space in the matrix space.
For the impatient reader, we point out that angles between these spaces are related to the behavior of
the staircase algorithm; note that R is always orthogonal to S. (We use ! \Delta; \Delta ? to represent the angle
between two spaces.)
angles components
A Staircase fails !
no weak stair no large large =2 small small
stair no large small =2 small large
small =2 large large
Here, by a weak stair [16], we mean the near rank deficiency of any superdiagonal block of the strictly
block upper triangular matrix A.
2. A Staircase Algorithm Failure to Motivate the Theory. Consider the two matrices
where ffi =1.5e-9 is approximately on the order of the square root of the double precision machine
roughly 2.2e-16. Both of these matrices clearly have the Jordan structure J 3 (0), but the
staircase algorithm on A 1 and A 2 can behave very differently.
To test this, we used the GUPTRI [11] algorithm. GUPTRI
1 requires an input matrix A and two
tolerance parameters EPSU and GAP. We ran GUPTRI on ~
2.2e-14 is roughly 100 times the double precision machine ffl. The singular values of each of the
two matrices ~
A 1 and ~
A 2 are oe 8.8816e-15. We set GAP to
be always  1, and let
we vary the value of a (The tolerance is effectively
a). Our observations are tabulated below.
a computed Jordan Structure for ~
A 1 computed Jordan Structure for ~
a
Here, we use J k () to represent a k \Theta k Jordan block with eigenvalue . In the table, typically
ff 6= fi 6= 0. Setting a small (smaller than here, which is the smaller singular value
in the second stage), the software returns two nonzero singular values in the first and second stages
of the algorithm and one nonzero singular value in the third stage. Setting EPSU \Theta GAP large (larger
than oe 2 here), we zero two singular values in the first stage and one in the second stage giving the
structure J 2 (0) \Phi J 1 (0) for both ~
A 1 and ~
(There is a matrix within O(10 \Gamma9 ) of A 1 and A 2 of the form
(0)). The most interesting case is in between. For appropriate EPSU \Theta GAP  a (between fl
and oe 2 here), we zero one singular value in each of the three stages, getting a J 3 (0) which is O(10 \Gamma14 )
away for A 2 , while we can only get a J 3 (0) which is O(10 \Gamma6 ) away for A 1 . In other words, the staircase
algorithm fails for A 1 but not for A 2 . As pictured in Figure 2.1, the A 1 example indicates that a matrix
1 GUPTRI [10, 11] is a "greedy" algorithm with a sophisticated thresholding procedure based on two input parameters
EPSU and GAP  1. We threshold oe (Defining oe n+1 j 0). The first
argument of the maximum oe k ensures a large gap between thresholded and non-thresholded singular values. The second
argument ensures that oe k\Gamma1 is small. Readers who look at the GUPTRI software should note that singular values are
ordered from smallest to largest, contrary to modern convention.
of the correct Jordan structure may be within the specified tolerance, but the staircase algorithm may
fail to find it.
Consider the situation when A 1 and A 2 are transformed using a random orthogonal matrix Q. As
a second experiment, we pick
\Gamma:39878 :20047 \Gamma:89487
\Gamma:84538 \Gamma:45853 :27400
and take ~
This will impose a perturbation of order ffl. We
ran GUPTRI on these two matrices; the following is the result:
a computed Jordan Structure for ~
A 1 computed Jordan Structure for ~
a
In the table, other values are the same as in the previous table.
In this case, GUPTRI is still able to detect a J 3 structure for ~
although the one it finds is O(10 \Gamma6 )
away. But it fails to find any J 3 structure at all for ~
A 1 . The comparison of A 1 and A 2 in the two
experiments indicates that the explanation is more subtle than the notion of a weak stair (a superdiagonal
block that is almost column rank deficient) [16].
In this paper we present a geometrical theory that clearly predicts the difference between A 1 and A 2 .
The theory is based on how close certain directions that we will denote staircase invariant directions
are to the tangent space of the manifold of matrices similar to the matrix with specified canonical form.
It turns out that for A 1 , these directions are nearly in the tangent space, but not for A 2 . This is the
crucial difference!
The tangent directions and the staircase invariant directions combine to form a "versal deformation"
in the sense of Arnold [1], but one with more useful properties for our purposes.
3. Staircase Invariant Space and Versal Deformations.
3.1. The Staircase Invariant Space and Related Subspaces. We consider block matrices
as in

Figure

3.1. Dividing a matrix A into blocks of row and column sizes we obtain a
general block matrix. A block matrix is conforming to A if it is also partitioned into blocks of
JA
Fig. 2.1. The staircase algorithm fails to find A1 at distance 2.2e-14 from ~
but does find a J3 (0) or a J2 (0) \Phi J1 (0)
if given a much larger tolerance. (The latter is ffi away from ~
.)
in the same manner as A. If a general block matrix has non-zero entries only in the
upper triangular blocks excluding the diagonal blocks, we call it a block strictly upper triangular
matrix. If a general block matrix has non-zero entries only in the lower triangular blocks including the
diagonal blocks, we call it a block lower triangular matrix. A matrix A is in staircase form if we
can divide A into blocks of sizes A is a strictly block upper triangular matrix
and every superdiagonal block has full column rank. If a general block matrix only has nonzero entries
on its diagonal blocks, and each diagonal block is an orthogonal matrix, we call it a block diagonal
orthogonal matrix. We call the matrix e B a block orthogonal matrix (conforming to
a block anti-symmetric matrix (conforming to (i.e. B is anti-symmetric with zero diagonal blocks.
Here, we abuse the word "conforming" since e B does not have a block structure.)
Definition 1. Suppose A is a matrix in staircase form. We call S a staircase invariant matrix
of A if S T is block lower triangular. We call the space of matrices consisting of all such S
the staircase invariant space of A, and denote it by S.
We remark that the columns of S will not be independent except possibly when can be
the zero matrix as an extreme case. However the generic sparsity structure of S may be determined by
general block matrix
block strictly upper
triangular matrix
block lower
triangular matrix00000000000000000000000000000011111111111111111111111111111111111111110000000000000011111111111111000000011111111111111
matrix in staircase
block diagonal
orthogonal matrix
block orthogonal
matrix
arbitrary block0000000000000000000011111111111111111111
full column rank block00000000000000111111111111111111111
orthogonal block special block zero block
Fig. 3.1. A schematic of the block matrices defined in the text.
the sizes of the blocks. For example, let A have the staircase form
\Theta \Theta
\Theta \Theta
\Theta \Theta
\Theta \Theta
\Theta \Theta
\Theta \Theta
\Theta
\Theta
\Theta
\Theta \Theta
\Theta \Theta
\Theta
\Theta
\Theta
\Theta \Theta \Theta
\Theta \Theta \Theta
\Theta \Theta \Theta
\Theta \Theta \Theta
\Theta \Theta \Theta
\Theta \Theta \Theta
\Theta \Theta \Theta
\Theta \Theta
\Theta \Theta
\Theta \Theta
\Theta \Theta
\Theta \Theta \Theta \Theta \Theta \Theta \Theta \ThetaC C C C C C C C A
is a staircase invariant matrix of A if every column of S is a left eigenvector of A. Here, the ffi notation
indicates 0 entries in the block lower triangular part of S that are a consequence of the requirement that
every column be a left eigenvector. This may be formulated as a general rule: if we find more than one
block of size n i \Theta n i then only those blocks on the lowest block row appear in the sparsity structure of
S. For example, the ffi do not appear because they are above another block of size 2. As a special case,
if A is strictly upper triangular, then S is 0 above the bottom row as is shown below. Readers familiar
with Arnold's normal form will notice that if A is a given single Jordan block in normal form, then S
contains the versal directions.
\Theta \Theta \Theta \Theta \Theta \Theta
\Theta \Theta \Theta \Theta \Theta
\Theta \Theta \Theta \Theta
\Theta \Theta \Theta
\Theta \Theta
\Theta \Theta \Theta \Theta \Theta \Theta \ThetaC C C C C A
Definition 2. Suppose A is a matrix. We call O(A) is a non-singular matrixg
the orbit of a matrix A. We call T any matrixg the tangent space of O(A) at
A.
Theorem 1. Let A be an n \Theta n matrix in staircase form, then the staircase invariant space S of A
and the tangent space T form an oblique decomposition of n \Theta n matrix space, i.e. R
Proof:
Assume that A i;j , the (i; j) block of A, is n i \Theta n j for and of course A
There are n 2
1 degrees of freedom in the first block column of S because there are n 1 columns and
each column may be chosen from the n 1 dimensional space of left eigenvectors of A. Indeed there are n 2
degrees of freedom in the ith block, because each of the n i columns may be chosen from the n i dimensional
space of left eigenvectors of the matrix obtained from A by deleting the first rows and columns.
The total number of degrees of freedom is
i , which combined with dim(T
gives the dimension of the whole space n 2 .
If S 2 S is also in T then S has the form AX \Gamma XA for some matrix X. Our first step will be
to show that X must have block upper triangular form after which we will conclude that AX \Gamma XA is
strictly block upper triangular. Since S is block lower triangular, it will then follow that if it is also in
must be 0.
Let i be the first block column of X which does not have block upper triangular structure. Clearly
the ith block column of XA is 0 below the diagonal block, so that the ith block column of
contains vectors in the column space of A. However every column of S is a left eigenvector of A from the
definition (notice that we do not require these column vectors of S to be independent, the one Jordan
block case is a good example.), and therefore orthogonal to the column space of A. Thus the ith block
column of S is 0, and from the full column rank conditions on the superdiagonal blocks of A, we conclude
that X is 0 below the block diagonal.
Definition 3. Suppose A is a matrix. We call O is a block anti-symmetric
matrix conforming to Ag the block orthogonal-orbit of a matrix A. We call
is a block anti-symmetric matrix conforming to Ag the block tangent space of the block
orthogonal orbit O b (A) at A. We call R j f block strictly upper triangular matrix conforming to Ag the
strictly upper block space of A.
Note that because of the complementary structure of the two matrices R and S, we can see that S
is always orthogonal to R.
Theorem 2. Let A be an n \Theta n matrix in staircase form, then the tangent space T of the orbit
O(A) can be split into the block tangent space T b of the orbit O b (A) and the strictly upper block space
Proof:
We know that the tangent space T of the orbit at A has dimension
. If we decompose
into a block upper triangular matrix and a block anti-symmetric matrix, we can decompose every
strictly upper triangular matrix and a matrix in T b . Since R, each of T b
and R has dimension  1=2(n
must both be exactly of dimension 1=2(n
Thus we know that they actually form a decomposition of T , and the strictly upper block space R can
also be represented as R j conforming to Ag:
Corollary 1. R n 2
Figure 3.2.
In Definition 3, we really do not need the whole set
merely need a small neighborhood around Readers may well wish to skip ahead to Section 4, but
for those interested in mathematical technicalities we review a few simple concepts. Suppose that we
have partitioned An orthogonal decomposition of n-dimensional space into k mutually
orthogonal subspaces of dimensions is a point on the flag manifold. (When this is
the Grassmann manifold). Equivalently, a point on the flag manifold is specified by a filtration, i.e.,
a nested sequence of subspaces V i of dimension
The corresponding decomposition can be written as
This may be expressed concretely. If from a unitary matrix U , we only define V i for
A
R
Fig. 3.2. A diagram of the orbits and related spaces. The similarity orbit at A is indicated by a surface O(A), the
block orthogonal orbit is indicated by a curve O b (A) on the surface, the tangent space of O b (A), T b is indicated by a line,
R which lies on O(A) is pictured as a line too, and the staircase invariant space S is represented by a line pointing away
from the plane.
the span of the first n 1 a point on the flag
manifold. Of course many unitary matrices U will correspond to the same flag manifold point. In an
open neighborhood of fe B g, near the point e the map between fe B g and an open subset of the
flag manifold is a one to one homeomorphism. The former set is referred to as a local cross section [24,
Lemma 4.1, page 123] in Lie algebra. No two unitary matrices in a local cross section would have the
same sequence of subspaces
3.2. Staircase as a Versal Deformation. Next, we are going to build up the theory of our versal
form. Following Arnold [1], a deformation of a matrix A, is a matrix A() with entries that are power
series in the complex variables  i , where convergent in a neighborhood of
with A.
A good introduction to versal deformations may be found in [1, Section 2.4] or [17]. The key property
of a versal deformation is that it has enough parameters so that no matter how the matrix is perturbed,
it may be made equivalent by analytic transformations to the versal deformation with some choice of
parameters. The advantage of this concept for a numerical analyst is that we might make a rounding
error in any direction and yet still think of this as a perturbation to a standard canonical form.
Let ae M be a smooth submanifold of a manifold M . We consider a smooth mapping A :   !M
of another manifold   into M , and let  be a point in   such that A() 2 N . The mapping A is called
transversal to N at  if the tangent space to M at A() is the direct sum
Here, TM A() is the tangent space of M at A(), TN A() is the tangent space of N at A(), T
is the tangent space of   at  and A   is the mapping from T    to TM A() induced by A (It is the
Jacobian).
Theorem 3. Suppose A is in staircase form. Fix S i 2
dim(S). It follows that
is a versal deformation of every particular A() for  small enough. A() is miniversal at
is a basis of S.
Proof:
Theorem 1 tells us the mapping A() is transversal to the orbit at A. From the equivalence of transversality
and versality [1], we know that A() is a versal deformation of A. Since the dimension of the
staircase invariant space S is the codimension of the orbit, A() given by Equation (3.1) is a miniversal
deformation if the S i are a basis for S (i.e. dim(S)). More is true, A() is a versal deformation
of every matrix in a neighborhood of A, in other words, the space S is transversal to the orbit of every
A(). Take a set of matrices X i s.t. the X i A \Gamma AX i form a basis of the tangent space T of the orbit at
A. We know T \Phi
, here \Phi implies T " so there is a fixed minimum angle ' between T
and S. For small enough , we can guarantee that the X i are still linearly independent
of each other and they span a subspace of the tangent space at A() that is at least, say, '=2 away from
S. This means that the tangent space at A() is transversal to S.
Arnold's theory concentrates on general similarity transformations. As we have seen above, the
staircase invariant directions are a perfect versal deformation. This idea can be refined to consider
similarity transformations that are block orthogonal. Everything is the same as above, except that we
add the block strictly upper triangular matrices R to compensate for the restriction to block orthogonal
matrices. We now spell this out in detail:
Definition 4. If the matrix C() is block orthogonal for every , then we refer to the deformation
as a block orthogonal deformation.
We say that two deformations A() and B() are block orthogonally-equivalent if there exists
a block orthogonal deformation C() of the identity matrix such that
We say that a deformation A() is block orthogonally-versal if any other deformation B()
is block orthogonally-equivalent to the deformation A(OE()). Here, OE is a mapping analytic at 0 with
Theorem 4. A deformation A() of A is block orthogonally-versal iff the mapping A() is transversal
to the block orthogonal-orbit of A at
Proof:
The proof follows Arnold [1, Sections 2.3 and 2.4] except that we use the block orthogonal version of the
relevant notions, and we remember that the tangents to the block orthogonal group are the commutators
of A with the block anti-symmetric matrices.
Since we know that T can be decomposed into T b \Phi R, we get:
Theorem 5. Suppose a matrix A is in staircase form. Fix S i 2
and k  dim(S). Fix R It follows that
is a block orthogonally-versal deformation of every particular A() for  small enough. A() is block
orthogonally-miniversal at A if fS i g, fR j g are bases of S and R.
It is not hard to see that the theory we set up for matrices with all eigenvalues 0 can be generalized
to a matrix A with different eigenvalues. The staircase form is a block upper triangular matrix, each
of its diagonal blocks of the form  staircase form defined at the beginning of this
chapter, and superdiagonal blocks arbitrary matrices. Its staircase invariant space is spanned by the
block diagonal matrices, each diagonal block being in the staircase invariant space of the corresponding
diagonal block A i . R space is spanned by the block strictly upper triangular matrices s.t. every diagonal
block is in the R space of the corresponding A i . T b is defined exactly the same as in the one eigenvalue
case. All our theorems are still valid. When we give the definitions or apply the theorems, we do not
really use the values of the eigenvalues, all that is important is how many different eigenvalues A has.
In other words, we are working with bundle instead of orbit.
These forms are normal forms that have the same property as the Arnold's normal form: they are
continuous under perturbation. The reason that we introduce block orthogonal notation is that the
staircase algorithm is a realization to first order of the block orthogonally-versal deformation, as we will
see in the next section.
4. Application to Matrix Staircase Forms. We are ready to understand the staircase algorithm
described in Section 1.2. We concentrate on matrices with all eigenvalues 0, since otherwise, the staircase
algorithm will separate other structures and continue recursively.
We use the notation stair(A) to denote the output A of the staircase algorithm as described in
Section 1.2. Now suppose that we have a matrix A which is in staircase form. To zeroth order, any
instance of the staircase algorithm replaces A with "
diagonal orthogonal.
Of course this does not change the staircase structure of A; the Q 0 represents the arbitrary rotations
within the subspaces, and can depend on how the software is written, and the subtlety of roundoff errors
when many singular values are 0. Next, suppose that we perturb A by fflE. According to Corollary 1, we
can decompose the perturbation matrix uniquely as
Theorem 6 states that in addition to some block diagonal matrix Q 0 , the staircase algorithm will apply
a block orthogonal similarity transformation to kill the perturbation in
Theorem 6. Suppose that A is a matrix in staircase form and E is any perturbation matrix. The
staircase algorithm (without zeroing) on A + fflE will produce an orthogonal matrix Q (depending on ffl)
and the output matrix
A has the same
staircase structure as A, "
S is a staircase invariant matrix of "
A and "
R is a block strictly upper triangular
matrix. If singular values are zeroed out, then the algorithm further kills "
S and outputs "
R+ o(ffl).
Proof:
After the first stage of the staircase algorithm, the first block column is orthogonal to the other columns,
and this property is preserved through the completion of the algorithm. Generally, after the ith iteration,
the ith block column below (including) the diagonal block is orthogonal to all other columns to its right,
and this property is preserved all through. So when the algorithm terminates, we will have a matrix
whose columns below (including) the diagonal block are orthogonal to all the columns to the right, in
other words, it is a matrix in staircase form plus a staircase invariant matrix.
We can always write the similarity transformation matrix as
a block diagonal orthogonal matrix and X is a block anti-symmetric matrix that does not depend on ffl
because of the local cross section property that we mentioned at the beginning of Section 3. Notice that
is not a constant matrix decided by A, it depends on fflE to its first order, we should have written
instead of Q 0 . However, we do not expand Q 0 since as long as it is a block diagonal
orthogonal transformation, it does not change the staircase structure of the matrix. Hence, we get
R and "
T b are respectively Q T
. It is easy to check that
T b is still in the space of "
A. X is a block anti-symmetric matrix satisfying "
AX.
We know that X is uniquely determined because the dimensions of "
b and the block anti-symmetric
matrix space are the same. The reason that "
AX hence the last equality in (4.1) holds is
because the algorithm forces the output form as described in the first paragraph of this proof: "
R
is in staircase form and ffl "
S is a staircase invariant matrix. Since (S \Phi R) " T b is the zero matrix, the T b
term must vanish.
To understand more clearly what this observation tells us, let us check some simple situations. If
the matrix A is only perturbed in the direction S or R, then the similarity transformation will be simply
a block diagonal orthogonal matrix Q 0 . If we ignore this transformation which does not change any
structure, we can think of the output to be unchanged from the input, this is the reason we call S the
staircase invariant space. The reason we did not include R into the staircase invariant space is that
fflR is still within O b (A). If the matrix A is only perturbed along the block tangent direction T b ,
then the staircase algorithm will kill the perturbation and do a block diagonal orthogonal similarity
transformation.
Although the staircase algorithm decides this Q 0 step by step all through the algorithm (due to SVD
rank decisions), we can actually think of the Q 0 as decided at the first step. We can even ignore this Q 0
because the only reason it comes up is that the svd we use follows a specific way to sort singular values
when they are different, and to choose the basis of the singular vector space when the same singular
values appear.
We know that every matrix A can be reduced to a staircase form under an orthogonal transformation,
in other words, we can always think of any general matrix M as P T AP , where A is in staircase form.
Thus in general, the staircase algorithm always introduces an orthogonal transformation and returns a
matrix in staircase form and a first order perturbation in its staircase invariant direction, i.e. stair(M
It is now obvious that if a staircase form matrix A has its S and T almost normal to each other,
then the staircase algorithm will behave very well. On the other hand, if S is very close to T then it
will fail. To emphasize this, we write it as a conclusion.
Conclusion 1. The angle between the staircase invariant space S and the tangent space T decides
the behavior of the staircase algorithm. The smaller the angle, the worse the algorithm behaves.
In the one Jordan block case, we have an if-and-only-if condition for S to be near T .
Theorem 7. Let A be an n \Theta n matrix in staircase form and suppose that all of its block sizes are
1 \Theta 1, then S(A) is close to T (A) iff the following two conditions hold:
(1)(row condition) there exists a non-zero row in A s.t. every entry on this row is o(1);
(2)(chain condition) there exists a chain of length with the chain value O(1), where k is the lowest
row satisfying (1).
Here, we call A i 1 ;i 2
a chain of length t and the product A i 1 ;i 2
is the
chain value.
Proof
Notice that S being close to T is equivalent to S being almost perpendicular to N , the normal space
of A. In this case, N is spanned by fI; A T ; A consists of matrices with nonzero
entries only in the last row. Considering the angle between any two matrices from the two spaces, it is
straightforward to show that S is almost perpendicular to N is equivalent to
(1) there exists a k s.t. the (n; entry of each of the matrices I; A
(2) if the entry is o(1), then it must have some other O(1) entry in the same matrix. Assume k is
the largest choice if there are different k's. By a combinatorial argument, we can show that these two
conditions are equivalent to the row and chain conditions respectively in our theorem.
Remark 1. Note that there exists an O(1) entry in a matrix is equivalent to say that there exists
a singular value of the matrix of O(1). So, the chain condition is the same as saying that the singular
values of A n\Gammak are not all O(ffl) or smaller.
Generally, we do not have an if-and-only-if condition for S to be close to T , we only have a necessary
condition, that is, only if at least one of the superdiagonal blocks of the original unperturbed matrix has
a singular value almost 0, i.e. it has a weak stair, will S be close to T . Actually, it is not hard to show
that the angle between T b and R is at most in the same order as the smallest singular value of the weak
stair. So, when the perturbation matrix E is decomposed into R are typically very
large, but whether S is large or not depends on whether S is close to T or not.
Notice that equation (4.1) is valid for sufficiently small ffl. What range of ffl is "sufficiently small"?
Clearly, ffl has to be smaller than the smallest singular value ffi of the weak stairs. Moreover, the algorithm
requires the perturbation along T and S to be both smaller than ffi. Assume the angle between T and
S is ', then generally, when ' is large, we would expect an ffl smaller than ffi to be sufficiently small.
However, when ' is close to 0, for a random perturbation, we would expect an ffl in the order of ffi=' to be
sufficiently small. Here, again, we can see that the angle between S and T decides the range of effective
ffl. For small ', when ffl is not sufficiently small, we observed some discontinuity in the 0th order term in
Equation (4.1) caused by the ordering of singular values during certain stages of the algorithm. Thus,
instead of the identity matrix, we get a permutation matrix in the 0th order term.
The theory explains why the staircase algorithm behaves so differently on the two matrices A 1 and
A 2 in Section 2. Using Theorem 7, we can see that A 1 is a staircase failure 2 is not
1). By a direct calculation, we find that the tangent space and the staircase invariant space of A 1 is
very close
this is not the situation for A 2
3).
When transforming to get ~
A 1 and ~
A 2 with Q, which is an approximate orthogonal matrix up to the
order of square root of machine precision ffl m , another error in the order of
it is comparable with ffi in our experiment, so the staircase algorithm actually runs on a shifted version
That is why we see R as large as an O(10 \Gamma6 ) added to J 3 in the second table
for ~
A 2 . We might as well call A 2 a staircase failure in this situation, but A 1 suffers a much worse failure
under the same situation, in that the staircase algorithm fails to detect a J 3 structure at all. This is
because the tangent space and the staircase invariant space are so close that the S and T component
are very large hence Equation (4.1) does not apply any more.
5. A Staircase Algorithm Failure to Motivate the Theory for Pencils. The pencil analog
to the staircase failure in Section 2 is
1.5e-8. This is a pencil with the structure L 1 \Phi J 2 (0). After we add a random perturbation
of size 1e-14 to this pencil, GUPTRI fails to return back the original pencil no matter which EPSU we
choose. Instead, it returns back a more generic L 2 \Phi J 1 (0) pencil O(ffl) away.
On the other hand, for another pencil with the same L 1 \Phi J 2 (0) structure:
GUPTRI returns an L 1 \Phi J 2 (0) pencil O(ffl) away.
At this point, readers may correctly expect that the reason behind this is again the angle between
two certain spaces as in the matrix case.
6. Matrix Pencils. Parallel to the matrix case, we can set up a similar theory for the pencil
case. For simplicity, we concentrate on the case when a pencil only has L-blocks and J(0)-blocks.
Pencils containing L T -blocks and non-zero (including 1) eigenvalue blocks can always be reduced to
the previous case by transposing and exchanging the two matrices of the pencil and/or shifting.
6.1. The Staircase Invariant Space and Related Subspaces for Pencils. A pencil (A; B) is
in staircase form if we can divide both A and B into block rows of sizes r columns
of sizes s strictly block upper triangular with every superdiagonal block having full
column rank and B is block upper triangular with every diagonal block having full row rank and the
rows orthogonal to each other. Here we allow s k+1 to be zero. A pencil is called conforming to (A; B)
if it has the same block structure as (A; B). A square matrix is called row (column) conforming to
if it has diagonal block sizes the same as the row (column) sizes of (A; B).
Definition 5. Suppose (A; B) is a pencil in staircase form and B d is the block diagonal part of
B. We call (SA ; SB ) a staircase invariant pencil of (A; B) if S T
has complimentary structure to (A; B). We call the space consisting of all such (SA ; SB ) the staircase
invariant space of (A; B), and denote it by S.
For example, let (A; B) have the staircase form
\Theta \Theta
\Theta \Theta
\Theta \Theta
\Theta \Theta
\Theta
\Theta
\Theta \Theta
\Theta \Theta
\Theta
\Theta
\Theta
\Theta \Theta \Theta
\Theta \Theta \Theta
\Theta \Theta
\Theta \Theta
\Theta \Theta
\Theta \Theta
\Theta
\Theta
\Theta \Theta
\Theta \Theta
\Theta \Theta
\Theta \Theta
\Theta
\Theta
\Theta \Theta
\Theta \Theta
\Theta
\Theta
then
\Theta \Theta \Theta
\Theta \Theta \Theta
\Theta \Theta
\Theta \Theta
\Theta \Theta
\Theta \Theta
\Theta \Theta \Theta \Theta \Theta \Theta \Theta \Theta7 7 7 7 7 7 7 5
\Theta \Theta \Theta
\Theta \Theta \Theta
\Theta \Theta \Theta
\Theta \Theta \Theta
\Theta \Theta \Theta ffi ffi
is a staircase invariant pencil of (A; B) if every column of SA is in the left null space of A and every
row of SB is in the right null space of B. Notice that the sparsity structure of SA and SB is at most
complimentary to that of A and B respectively, but SA and SB are often less sparse, because of the
requirement on the nullspace. To be precise, if we find more than one diagonal block with the same
size, then among the blocks of this size, only the blocks on the lowest block row appear in the sparsity
structure of SA . If any of the diagonal blocks of B is a square block, then SB has all zero entries
throughout the corresponding block column.
As special cases, if A is a strictly upper triangular square matrix and B is an upper triangular square
matrix with diagonal entries nonzero, then SA only has nonzero entries in the bottom row and SB is
simply a zero matrix. If A is a strictly upper triangular n \Theta (n + 1) matrix and B is an upper triangular
with diagonal entries nonzero, then (SA ; SB ) is the zero pencil.
Definition 6. Suppose (A; B) is a pencil. We call O(A; B) are non-singular
square matrices g the orbit of a pencil (A; B). We call T are any square
matrices g the tangent space of O(A; B) at (A; B).
Theorem 8. Let (A; B) be an m \Theta n pencil in staircase form, then the staircase invariant space S of
B) and the tangent space T form an oblique decomposition of m \Theta n pencil space, i.e. R
Proof:
The proof of the theorem is similar to that of Theorem 1; first we prove the dimension of S(A; B) is the
same as the codimension of T (A; B), then we prove induction. The readers may try to
fill out the details.
Definition 7. Suppose (A; B) is a pencil. We call O b (A; B) j fP is a block
anti-symmetric matrix row conforming to (A; B); is a block anti-symmetric matrix column
conforming to (A; B) the block orthogonal-orbit of a pencil (A; B). We call T b j
X is a block anti-symmetric matrix row conforming to (A; B), Y is a block anti-symmetric matrix column
conforming to (A; B)g the block tangent space of the block orthogonal-orbit O b (A; B) at (A; B). We
call R j U is a block upper triangular matrix row conforming to (A; B), V is a
block upper triangular matrix column conforming to (A; B)g the block upper pencil space of (A; B).
Theorem 9. Let (A; B) be an m \Theta n pencil in staircase form, then the tangent space T of the orbit
O(A; B) can be split into the block tangent space T b of the orbit O b (A; B) and the block upper pencil
space R, i.e.
Proof:
This can be proved by a very similar argument concerning the dimensions as for matrix, in which the
dimension of R is 2
the dimension of T b is
the codimension
of the orbit O(A; B) (or T ) is
Corollary 2. R
6.2. Staircase as a Versal Deformation for pencils. The theory of versal forms for pencils
[17] is similar to the one for matrices. A deformation of a pencil (A; B) is a pencil (A; B)() with
entries power series in the real variables  i . We say that two deformations (A; B)() and (C; D)() are
equivalent if there exist two deformations P () and Q() of identity matrices such that (A;
Theorem 10. Suppose (A; B) is in staircase form. Fix S i 2
dim(S). It follows that
is a versal deformation of every particular (A; B)() for  small enough. (A; B)() is miniversal at
is a basis of S.
Definition 8. We say two deformations (A; B)() and (C; D)() are block orthogonally-
equivalent if there exist two block orthogonal deformations P () and Q() of the identity matrix such
that are exponentials of matrices which are
conforming to (A; B) in row and column respectively.
We say that a deformation (A; B)() is block orthogonally-versal if any other deformation
(C; D)() is block orthogonally-equivalent to the deformation (A; B)(OE()). Here, OE is a mapping holomorphic
at 0 with
Theorem 11. A deformation (A; B)() of (A; B) is block orthogonally-versal iff the mapping
(A; B)() is transversal to the block orthogonal-orbit of (A; B) at
This is the corresponding result to Theorem 4.
Since we know that T can be decomposed into T b \Phi R, we get:
Theorem 12. Suppose a pencil (A; B) is in staircase form. Fix S i 2
S and k  dim(S). Fix R It follows that
is a block orthogonally-versal deformation of every particular (A; B)() for  small enough. (A; B)()
is block orthogonally-miniversal at (A; B) if fS i g, fR j g are bases of S and R.
Notice that as in the matrix case, we can also extend our definitions and theorems to the general form
containing L T -blocks and non-zero eigenvalue blocks, and again, we will not specify what eigenvalues
they are and hence get into the bundle case. We only want to point out one particular example here.
If (A; B) is in the staircase form of will be a strictly upper triangular matrix with
nonzero entries on the super diagonal and B will be a triangular matrix with nonzero entries on the
diagonal except the (n will be the zero matrix and SB will be a matrix with the
only nonzero entry on its (n
7. Application to Pencil Staircase Forms. We concentrate on L \Phi J(0) structures only, since
otherwise, the staircase algorithm will separate all other structures and continue similarly after a shift
and/or transpose on that part only. As in the matrix case, the staircase algorithm basically decomposes
the perturbation pencil into three spaces T b , R, and S and kills the perturbation in T b .
Theorem 13. Suppose that (A; B) is a pencil in staircase form and E is any perturbation pencil.
The staircase algorithm (without zeroing) on (A; B)+ fflE will produce two orthogonal matrices P and Q
(depending on ffl) and the output pencil stair((A; B)+
R)+o(ffl),
B) has the sane staircase structure as
S is a staircase invariant pencil of ( "
B) and
R is in the block upper pencil space R. If singular values are zeroed out, then the algorithm further kills
S and output
We use a formula to explain the statement more clearly:
Similarly, we can see that when a pencil has its T and S almost normal to each other, the staircase
algorithm will behave well. On the other hand, if S is very close to T , then it will behave badly. This
is exactly the situation in the two pencil examples in Section 5. Although the two pencils are both
ill conditioned, a direct calculation shows that the first pencil has its staircase invariant space very
close to the tangent space (the angle !
while the second one does not (the angle
The if-and-only-if condition for S to be close to T is more difficult than in the matrix case. One
necessary condition is that one super diagonal block of A is almost of not full column rank or one diagonal
block of B is almost not full row rank. This is usually referred to as weak coupling.
8. Examples: The geometry of the Boley pencil and others. Boley [3, Example 2, Page
639] presents an example of a 7 \Theta 8 pencil (A; B) that is controllable (has generic Kronecker structure)
yet it is known that an uncontrollable system (non-generic Kronecker structure) is nearby at a distance
6e-4. What makes the example interesting is that the staircase algorithm fails to find this nearby
uncontrollable system while other methods succeed. Our theory provides a geometrical understanding
of why this famous example leads to staircase failure: the staircase invariant space is very close to the
tangent space.
The pencil that we refer to is (A; B(ffl)), where
and
(The dots refer to zeros, and in the original Boley example
the staircase algorithm predicts a distance of 1, and is therefore off by nearly four
orders of magnitude. To understand the failure, our theory works best for smaller values of ffl, but it is
still clear that even for there will continue to be difficulties.
It is useful to express the pencil (A; B(ffl)) as is zero except for
a "one" in the (7,7) entry of its B part. P 0 is in the bundle of pencils whose Kronecker form is L 6 +J 1 (\Delta)
and the perturbation E is exactly in the unique staircase invariant direction (hence the notation "S")
as we pointed out at the end of Section 6.
The relevant quantity is then the angle between the staircase invariant space and the pencil space.
An easy calculation reveals that the angle is very small: ' radians. In order to get a feeling
for what range of ffl first order theory applies, we calculated the exact distance d(ffl) j d(P (ffl); bundle)
using the nonlinear eigenvalue template software [33]. To first order, Figure 8.1 plots the
distances first for ffl 2 [0; 2] and then a closeup for
size of perturbation
the
distance
to
the
orbit
size of perturbation
the
distance
to
the
orbit
Fig. 8.1. The picture to explain the change of the distance of the pencils P0 + fflE to the bundle of L6
changes. The second subplot is part of the first one at the points near
Our observation based on this data suggests that first order theory is good to two decimal places
for one place for ffl  10 \Gamma2 . To understand the geometry of staircase algorithmic failure,
one decimal place or even merely an order of magnitude is quite sufficient.
In summary, we see clearly that the staircase invariant direction is at a small angle to the tangent
space, and therefore the staircase algorithm will have difficulty finding the nearest pencil on the bundle
or predicting the distance. This difficulty is quantified by the angle ' S .
Since the Boley example is for we computed the distance well past 1. The breakdown of
first order theory is attributed to the curving of the bundle towards S. A three dimensional schematic
is portrayed in Figure 8.2.
Fig. 8.2. The staircase algorithm on the Boley example. The surface represents the orbit O(P0 ). Its tangent space
at the pencil P0 , T (P0 ), is represented by the plane on the bottom. P1 lies on the staircase invariant space S inside the
"bowl". The hyperplane of uncontrollable pencils is represented by the plane cutting through the surface along the curve
C. It intersects T (P0 ) along L. The angle between L and S is ' c . The angle between S and T (P0 ), ' S , is represented by
the angle "HP0P1 .
The relevant picture for control theory is a planar intersection of the above picture. In control theory,
we set the special requirement that the "A" matrix has the form [0 I]. Pencils on the intersection of
this hyperplane and the bundle are termed "uncontrollable."
We analytically calculated the angle ' c between S and the tangent space for the "uncontrollable
We found that ' 0:0040. Using the nonlinear eigenvalue template software [33], we
numerically computed the true distance from fflE to the "uncontrollable surfaces" and calculated
the ratio of this distance to ffl, we found that for ffl ! 8e \Gamma 4, the ratio agrees with ' very well.
We did a similar analysis on the three pencils C 1 , C 2 C 3 given by J. Demmel and B. Kagstrom [12].
We found that the sin values of the angles between S and T are respectively 2.4325e-02, 3.4198e-02
and 8.8139e-03, and the sin values between T b and R are respectively 1.7957e-02 7.3751e-03 and
3.3320e-06. This explains why we saw the staircase algorithm behave progressively worse on them.
Especially, it explains why when a perturbation about 10 \Gamma3 is added to these pencils, C 3 behaves
dramatically worse then C 1 and C 2 . The component in S is almost of the same order as the entries of
the original pencil.
So we conclude that the reason the staircase algorithm does not work well on this example is because
actually a staircase failure, in that its tangent space is very close to its staircase
invariant space and also the perturbation is so large that even if we know the angle in advance we can
not estimate the distance well.

Acknowledgement

. The authors thank Bo Kagstrom and Erik Elmroth for their helpful discussion
and their conlab software for easy interactive numerical testing. The staircase invariant directions were
originally discovered for single Jordan blocks with Erik Elmroth while he was visiting MIT during the
fall of 1996.



--R

On matrices depending on parameters.
An improved algorithm for the computation of Kronecker's canonical form of a singular pencil.
Estimating the sensitivity of the algebraic structure of pencils with simple eigenvalue estimates.
The algebraic structure of pencils and block Toeplitz matrices.
Placing zeroes and the Kronecker canonical form.
Lectures on Finite Precision Computations.
The dimension of matrices (matrix pencils) with given Jordan (Kronecker) canonical forms.
Stably computing the Kronecker structure and reducing subspace of singular pencils A
Computing stable eigendecompositions of matrix pencils.
The generalized Schur decomposition of an arbitrary pencil A
The generalized Schur decomposition of an arbitrary pencil A
Accurate solutions of ill-posed problems in control theory
The computation of Kronecker's canonical form of a singular pencil.
The generalized eigenstructure problem in linear system theory.
Reducing subspaces: definitions
Oral communication.
A geometric approach to perturbation theory of matrices and matrix pencils: Part 1: versal deformations.
A geometric approach to perturbation theory of matrices and matrix pencils: Part 2: stratification-enhanced staircase algorithm
The set of 2-by-3 matrix pencils - Kronecker structures and their transitions under perturbations
Computation of zeros of linear multivariable systems.
The application of singularity theory to the computation of Jordan Canonical Form.
Matrix Computations.

Differential Geometry
The generalized singular value decomposition and the general A

ALGORITHM 560: JNF
An algorithm for numerical computation of the Jordan normal form of a complex matrix.
Matrix Pencils
Robust pole assignment in linear state feedback.
On a method of solving the complete eigenvalue problem of a degenerate matrix.

Nonlinear eigenvalue problems.
An algorithm for numerical determination of the structure of a general matrix.
Computing the distance to an uncontrollable system.
--TR

--CTR
Naren Ramakrishnan , Chris Bailey-Kellogg, Sampling Strategies for Mining in Data-Scarce Domains, Computing in Science and Engineering, v.4 n.4, p.31-43, July 2002
kronecker structure;jordan structure;SVD;staircase algorithm;versal deformation