Matrix Powers in Finite Precision Arithmetic.
If $A$ is a square matrix with spectral radius less than 1 then $A^k \to 0$ as $k \to \infty$, but the powers computed in finite precision arithmetic may or may not converge. We derive a sufficient condition for $fl(A^k) \to 0$ as $k\to\infty$ and a bound on $\norm{fl(A^k)}$, both expressed in terms of the Jordan canonical form of $A$. Examples show that the results can be sharp. We show that the sufficient condition can be rephrased in terms of a pseudospectrum of $A$ when $A$ is diagonalizable, under certain assumptions. Our analysis leads to the rule of thumb that convergence or divergence of the computed powers of $A$ can be expected according as the spectral radius computed by any backward stable algorithm is less than or greater than 1.
Introduction
. Many numerical processes depend for their success upon the
powers of a matrix tending to zero. A fundamental example is stationary iteration
for solving a linear system b, in which a sequence of vectors is defined by
and M is nonsingular. The errors e
so the iteration converges for all x 0 if (M
Many theorems are available about the convergence of stationary iteration,
but virtually all of them are concerned with exact arithmetic (for exceptions see [12],
[13] and the references therein). While the errors in stationary iteration are not
precisely modelled by the errors in matrix powering, as matrix powers are not formed
explicitly, the behaviour of the computed powers f l((M can be expected to
give some insight into the behaviour of stationary iteration (indeed, the basic error
recurrences in [12] and [13] involve powers of M \Gamma1 N acting on vectors of rounding
errors).
In [18, Chapter 20], Ostrowski proves a theorem about a product of perturbed
states "assures the theoretical stability of the convergence
of A - to 0 with respect to rounding off" as - ! 1 for any matrix A with spectral
radius ae(A) ! 1. Although Ostrowski's theorem is correct, its interpretation with
respect to computed powers is not as simple as this statement implies, because for
any finite precision arithmetic, no matter how accurate, there are matrices that are
sensitive enough to perturbations to cause the theoretically convergent sequence of
powers to diverge. To illustrate this point, Figure 1.1 1 plots the 2-norms of the first
200 powers of a 14 \Theta 14 nilpotent matrix C 14 discussed by Trefethen and Trummer
[23] (see x3 for details). The plot confirms the statement of these authors that the
matrix is not power-bounded in floating point arithmetic, even though its 14th power
should be zero. The powers for our plot were computed in Matlab, which has unit
roundoff Reichel and Trefethen [19] also give an example
Revised version dated March 7, 1994. To appear in SIAM J. Matrix Anal. Appl., April 1995.
y Department of Mathematics, University of Manchester, Manchester, M13 9PL, England
(na.nhigham@na-net.ornl.gov). The work of this author was supported by Science and Engineering
Research Council grant GR/H52139.
z Department of Mathematics, University of Manchester, Manchester, M13 9PL, England.
Present address: Department of Mathematics, University of Strathclyde, Glasgow, G1 1XH, Scotland
(p.a.knight@strath.ac.uk). This author was supported by an SERC Research Studentship.
1 As in all our plots of norms of powers, k on the x-axis is plotted against kf l(A k )k 2 on the y-axis.
Fig. 1.1. Diverging powers of a nilpotent matrix, C 14 .
of a matrix that is nilpotent in theory but not power-bounded in practice. In this
paper we determine conditions on a matrix A that ensure that the computed powers
converge to zero.
In Section 2 we examine the behaviour of matrix powers in exact arithmetic. In
particular, we review a number of bounds on the norms of powers. In Section 3 we
use the Jordan canonical form of A to bound kf l(A k )k and to determine a sufficient
condition for 1. We also show that for certain matrices our
bounds are tight. Finally, in Section 4 we rephrase our sufficient condition in terms of
a pseudospectrum of A, under certain assumptions, including that A is diagonalizable;
the modified result is not any sharper than the original, but offers an alternative
viewpoint that is intuitively attractive.
In our analysis we use the standard model for floating point arithmetic:
where u is the unit roundoff. This model is valid for machines that do not use a guard
digit in addition and subtraction.
We will use the Frobenius norm,
and the p-norms
From x3 onwards
we will drop the subscripts on k \Delta k p , since all the norms from that point on are p-norms.
2. Matrix Powers in Exact Arithmetic. We begin by discussing the behaviour
of matrix powers in the absence of rounding errors. In exact arithmetic the
limiting behaviour of the powers of A 2 C n\Thetan is determined by A's eigenvalues. If the
spectral radius
if A has a defective eigenvalue - such that does
not converge if A has a nondefective eigenvalue - 6= 1 such that (although the
norms of the powers may converge); otherwise, the only eigenvalue of modulus 1 is the
nondefective eigenvalue 1, and A k converges to a nonzero matrix. These statements
are easily proved using the Jordan canonical form
where X is nonsingular and
We will call a matrix for which A k ! 0 as k !1 (or
equivalently, convergent matrix.
The norm of a convergent matrix can be arbitrarily large, as is shown trivially by
the example
While the spectral radius determines the asymptotic rate of growth of matrix pow-
ers, the norm influences the initial behaviour of the powers. The interesting result
that norm (see [14, p. 299], for example)
confirms the asymptotic role of the spectral radius. An important quantity is the
which can be arbitrarily large for a convergent matrix,
as can be seen from A 3 (ff), the 3 \Theta 3 analogue of the matrix in (2.2), for which
O(ff). Figure 2.1 shows an example of the hump phenomenon:
the plot is for A 3 (2) with 3:57. The shape of the plot is
typical of that for a convergent matrix with norm bigger than 1. Note that if A is
normal (so that in (2.1a) J is diagonal and X can be taken to be unitary) we have
so the problem of bounding kA k k is of interest
only for nonnormal matrices. The hump phenomenon arises in various areas of
numerical analysis. For example, it is discussed for matrix powers in the context of
stiff differential equations by D. J. Higham and Trefethen [8], and by Moler and Van
Loan [17] for the matrix exponential e At with t !1.
In the rest of this section we briefly survey bounds for kA k k. First, however, we
comment on the condition number that appears in various bounds
in this paper. The matrix X in the Jordan form (2.1a) is by no means unique [3,
pp. 220-221], [6]: if A has distinct eigenvalues (hence J is diagonal) then X can be
replaced by XD, for any nonsingular diagonal D, while if A has repeated eigenvalues
then X can be replaced by XT , where T is a block matrix with block structure
conformal with that of J and which contains some arbitrary upper trapezoidal Toeplitz
blocks. We adopt the convention that -(X) denotes the minimum possible value of
-(X) over all possible choices of X. In general it is difficult to determine this optimal
value. However, for any nonsingular X we have the bound
equality if there is a nonzero
ff such that Theorem 4.3.5]. If A has distinct eigenvalues
then this lower bound is the same for all X in the Jordan form and is the minimum
Fig. 2.1. A typical hump for a convergent, nonnormal matrix.
value of -F (X). An alternative approach for matrices with distinct eigenvalues is
to insist that the columns of X have unit 2-norm, for this gives a 2-norm condition
number within a factor n 1=2 of the minimum, in view of a result by van der Sluis on
diagonal scalings [24, Theorem 3.5]. However, we will see in x3 that to appreciate
fully the various instability phenomena, we must consider defective problems.
If A is diagonalizable then, from (2.1a), we have the bound
for any p-norm. (Since ae(A) - kAk for any norm, we also have the lower bound
.) This bound is unsatisfactory for two reasons. First, by choosing
A to have well-conditioned large eigenvalues and ill-conditioned small eigenvalues
we can make the bound arbitrarily pessimistic. Second, it models norms of powers
of convergent matrices as monotonically decreasing sequences, which is qualitatively
incorrect if there is a large hump.
The Jordan canonical form can also be used to bound the norms of the powers of
a defective matrix. If XJX \Gamma1 is the Jordan canonical form of ffi \Gamma1 A then
for all ffi ? 0. This is a special case of a result of Ostrowski [18, Theorem 20.1],
and a proof is straightforward: We can
is the off-diagonal part of the Jordan form. Then
taking norms. An alternative way of writing this bound
is
1). Note that this is not the same
X as in (2.4): multiplying A by a scalar changes -(X) when A is not diagonalizable.
Both bounds suffer from the same problems as the bound (2.3) for diagonalizable
matrices.
Another bound in terms of the Jordan canonical form (2.1) of A is given by
Gautschi [4]. For convergent matrices, it can be written in the form
and c is a constant depending only on A (c is not
defined explicitly in [4]). The factor k p\Gamma1 makes this bound somewhat more effective
at predicting the shapes of the actual curve than (2.4), but again c can be unsuitably
large.
Another way to estimate kA k k is to introduce a measure of nonnormality. Consider
the Schur decomposition Q H strictly upper triangular,
and let S represent the set of all such N . The nonnormality of A can be measured by
Henrici's departure from normality [7]
For the Frobenius norm, Henrici shows that kNkF is independent of the particular
Schur form and that
L'aszl'o [15] has recently shown that \Delta F (A) is within a constant factor of the distance
from A to the nearest normal matrix:
normalg. Henrici uses the departure from normality
to derive the 2-norm bounds
Empirical evidence suggests that the first bound in (2.6) can be very pessimistic.
However, for normal matrices both the bounds are equalities. A bound of the same
form as the first bound in (2.6), but with replacing and with an extra
is obtained from a bound of Stafney in [20, Theorem 2.1] for kp(A)k,
where p is a polynomial.
Another bound involving nonnormality is given by Golub and Van Loan [5,
Lemma 7.3.2]. They show that, in the above notation,
for any ' - 0. This bound is an analogue of (2.4) with the Schur form replacing the
Jordan form. Again, there is equality when A is normal (if we set
To compare bounds based on the Schur form with ones based on the Jordan form
we need to compare \Delta(A) with -(X). If A is diagonalizable then [16, Theorem 4]
F
and it can be shown by a 2 \Theta 2 example that minX - 2 (X) can exceed \Delta F (A)=kAkF
by an arbitrary factor [2, Section 8.1.2], [1, Section 4.2.7].
Another tool that can be used to bound the norms of powers is the pseudospec-
trum of a matrix [22]. The ffl-pseudospectrum of A 2 C n\Thetan is defined for a given ffl ? 0
to be the set
z is an eigenvalue of
and it can also be represented, in terms of the resolvent (zI
As Trefethen notes [22], by using the Cauchy integral representation of A k (which
involves a contour integral of the resolvent) one can show that
where the ffl-pseudospectral radius
This bound is very similar in flavour to (2.4). The difficulty is transferred from
estimating -(X) to choosing ffl and estimating ae ffl (A).
Finally, we mention that the Kreiss matrix theorem provides a good estimate of
sup k-0 kA k k for a general A 2 C n\Thetan , albeit in terms of an expression that involves
the resolvent and is not easy to compute:
references are given by Wegert and Trefethen [25].
3. Bounds for Finite Precision Arithmetic. The formulae A \Delta A k or A k \Delta A
can be implemented in several ways, corresponding to different loop orderings in
each individual product, but as long as each product is formed using the standard
variations satisfy the same rounding error
bounds. We do not analyse here the use of the binary powering technique, where,
for example, A 9 is formed as alternate multiplication on the left and
right: or the use of fast matrix multiplication techniques
such as Strassen's method, since none of these methods is equivalent to repeated
multiplication in finite precision arithmetic.
We suppose, without loss of generality, that the columns of A m are computed one
at a time, the jth as f is the jth unit vector. Standard
error analysis shows that the jth computed column of A m satisfies
where
with c n a constant of order n. (The inequality and absolute value are taken compo-
nentwise.) This bound holds for both real and complex matrices. It follows that
and so a sufficient condition for convergence of the computed powers is that
This result is useful in certain special cases: is triangular or has
a checkerboard sign pattern (since then
normal then ae(jAj) -
(this bound being attained for a Hadamard matrix);
and in Markov processes, where the a ij are transition probabilities, A. However,
in general ae(jAj) can exceed ae(A) by an arbitrary factor.
To obtain sharper and more informative results it is necessary to use more information
about the matrix. Although the Jordan form is usually avoided by numerical
analysts because of its sensitivity to perturbations, it is convenient to work with in
this application and leads to informative results.
We point out that, because the analysis below is based on (3.1), our proofs of
sufficient conditions for only trivial changes, sufficient conditions
for 0, for any vector b. These conditions do not, however, exploit
any special relations between A and b (such, as for example, b being an eigenvector
of A).
3.1. Nilpotent Matrices. We begin by considering nilpotent matrices, that is,
those whose spectral radius is zero. The fact that nth power of an n \Theta n nilpotent
matrix is zero simplifies the analysis. The following theorem gives a bound on the
norm of a computed power, together with a condition for the limit of the powers to
be zero.
Theorem 3.1. Let A 2 C n\Thetan be a nilpotent matrix with the Jordan canonical
form (2.1). A sufficient condition for
for some p-norm, where u is the unit roundoff and dn is a modest constant that
depends only on n. Furthermore, if, for some k - 1 and ' ? 1,
Proof. Taking norms in (3.1) we have
Using the inequality
from [10], we have
Expanding this product and collecting together terms with the same number of \DeltaA i
factors we obtain the bound
From (3.2) we find, using (3.6) and an analogous result involving duality, that
since
every term in the first r \Gamma 1 summations in (3.7) contains a factor
all these terms disappear. Furthermore, in the remaining summations we need only
count terms in which all the exponents of J are less than t (again, the other terms
disappear). Overall, we have, using the fact that
Now suppose that (3.4) holds with
n , for some ' ? 1. Then, for r - k,
This gives the second part of the theorem. The first part follows immediately by
choosing ffl, with ffl an arbitrarily small positive number, and taking the limit
as k !1.
In practice we may have a computed matrix b
A - A that is not exactly nilpotent.
As long as
nukAk, we can absorb the error A \Gamma b
A into the terms \DeltaA i in
the proof, and so by applying the theorem to A we will obtain conclusions valid for
A.
To exhibit the sharpness of the bounds we give the following example, using the
Chebyshev spectral differentiation matrix Cn 2 IR n\Thetan described in [23]. The matrix
Fig. 3.1. Converging powers of the nilpotent matrix C 8 .
Cn arises from degree polynomial interpolation of n arbitrary data values at n
Chebyshev points, including a boundary condition at x 1. It is nilpotent and is
similar to a single Jordan block of dimension n. We generate Cn in Matlab using the
routine chebspec from the test collection of Higham [9], [11].

Figure

3.1 shows the 2-norms of the computed powers of C 8 and Figure 1.1 those
of C 14 . The powers of C 8 converge to zero, while the powers of C 14 diverge.
To check the sharpness of the bounds in Theorem 3.1 we need an estimate of the
condition number of the matrix Xn in the Jordan canonical form of Cn . We outline our
approach in an appendix. Our estimate for

Table

3.1 gives the order of the bound (3.5) for a number of powers, with
' chosen as large as possible so that (3.4) is satisfied (we take of the
actual value dn - n 5=2 for this example, to allow for the inevitable overestimation of
errors inherent in a strict rounding error bound of this type). The actual computed
order is given for comparison and clearly there is reasonable agreement. According
to (3.3), we require dn-(X)kAku ! 1 to guarantee that the computed powers of A
converge to zero. For C 14 we have u- 2 (X)kC 14 k 2 - 0:28 so, allowing for dn , (3.3)
correctly does not predict convergence of the computed powers.
To emphasize that the behaviour of the computed powers is scale-dependent, we
mention that the computed powers of 15C 8 diverge. Again, this is in accord with
Theorem 3.1 because u- 2 (X)k15C 8 k 2 - 2:7. Finally, we note that for C 12 , Theorem
3.1 again correctly predicts convergence of the computed powers, but the powers
computed by alternate left and right multiplication and by binary powering diverge;
this confirms that our analysis is not applicable to these forms of multiplication.
3.2. General Matrices. Now we turn to general convergent matrices. In contrast
to the theory we have developed for nilpotent matrices, we need separate theorems
to describe the limiting behaviour of the matrix powers and to bound the norm
for a finite exponent. In the following theorem we give a sufficient condition, based
on the Jordan canonical form, for the computed powers of a matrix to converge to
zero.
Table
Expected and actual orders of kf l(C m)k 2 .
Power Bound Actual
Theorem 3.2. Let A 2 C n\Thetan with the Jordan form (2.1) have spectral radius
1. A sufficient condition for
for some p-norm, where is a modest constant depending only on
n.
Proof. Since any two p-norms differ by a factor at most n, we need only show
convergence for one particular norm. We choose the 1-norm.
It is easy to see that if we can find a nonsingular matrix S such that
for all i, then the product
In the rest of the proof we construct such a matrix
S for the \DeltaA i in (3.1).
Let
Now consider the matrix Its ith diagonal block is of the form -
where the only nonzeros in N are 1s on the first superdiagonal, and
so
Now we set and we determine ' so that (3.9) is
satisfied, that is, so that -1 From (3.2) and (3.10) we have
Therefore (3.9) is satisfied if
that is, if
Fig. 3.2. Diverging powers of
If the integer t is greater than 1 then the function has a maximum
on [0; 1] at '
We conclude that for all integers 1 - t - n,
is sufficient to ensure that (3.9) holds. Thus the theorem is proved, with
If A is normal then takes the
This condition is also easily derived by taking 2-norms in (3.1) and (3.2).
Again, we can show the sharpness of this bound by using the Chebyshev spectral
differentiation matrix Cn , this time adding multiples of the identity matrix.

Figure

3.2 shows the non-converging 2-norms of the first 200 computed powers of
0:36I. Since the same matrix X takes both C 13 and A to Jordan form,
we can use the same routines as for our nilpotent examples to estimate -(X). Our
estimate for - 2 13 in this case is 3.05. On the other hand, the
computed powers of converge to zero, and - 2
0:01. Thus our bound (3.8) is reasonably sharp.

Figure

3.2 reveals an interesting scalloping pattern in the curve of the norms.
In

Figures

1.1 and 3.1 for nilpotent matrices the norm dips whenever the power is
a multiple of the dimension of the matrix. Here the norm first dips for kf l(A 28 )k 2
and then regularly after every further 20 powers, but the point of first dip and the
dipping intervals can be altered by adding different multiples of the identity matrix.
The reason for this behaviour is not clear.
A difficulty we have when attempting to bound kf l(A m )k for a finite m is that,
as explained in Section 2, we do not have a good estimate of the true value
we do have such an estimate we can prove results similar to (3.4) and (3.5), although
we are not able to determine a precise bound for kf l(A m )k in simple form.
Theorem 3.3. Let A 2 C n\Thetan with the Jordan form (2.1) have spectral radius
q be such that kA q and suppose that, for some
Proof. We omit the proof of the theorem, which is very similar to the proof of
Theorem 3.1.
We conclude this section by commenting that the proof of Theorem 3.2 can be
adapted to use the Schur decomposition in place of the Jordan canonical form. The
modified analysis leads to the sufficient condition
for is the strictly upper triangular part of the Schur
form. This condition is weaker than (3.8) in two respects. First, it takes no account
of the defectiveness of A, because it contains a power n on the right-hand side instead
of under the scaling A / ffA the left-hand side of (3.11)
scales by jffj n , which tends to make the left-hand side of (3.11) much larger than
that of (3.8) when kAkF ? 1. It is an open question how to obtain a sharp sufficient
condition for convergence in terms of the Schur decomposition.
4. A Pseudospectral Approach. In this section we show how the pseudospec-
trum can be used to predict the limiting behaviour of a computed sequence of powers.

Figure

4.1 shows approximations to pseudospectra for the matrices in the examples
of

Figures

1.1, 3.1 and 3.2; we have approximated   ffl (A) with taking
simplicity. The inner ring is an approximation to the pseudospectrum of
C 8 , that of C 14 is marked by "+" and that of C 13 +0:36I is marked by "o". The solid
curve is the unit disc.
A heuristic argument based on (3.1) and (3.2) suggests that, if for randomly
chosen perturbations \DeltaA i with most of the eigenvalues of the
perturbed matrices lie outside the unit disc, then we can expect a high percentage of
the terms in (3.1) to have spectral radius bigger than one and hence we can
expect the product to diverge. On the other hand, if the c n ukAk-pseudospectrum is
wholly contained within the unit disc, each have spectral radius less than
one and the product can be expected to converge. (Note, however, that if ae(A)
and it is not necessarily the case that ae(AB) ! 1.) To make this heuristic
precise, we need an analogue of Theorem 3.2 phrased in terms of the pseudospectrum
rather than the Jordan form.
To obtain such an analogue directly from Theorem 3.2 we need to relate -(X) to
the pseudospectral radius ae ffl (A) (see (2.8)). If we can show that
for a particular ffl, then ae which is a condition
of the same form as (3.8). In Theorem 4.2 we show that (4.1) holds for diagonalizable
Fig. 4.1. Pseudospectra of the three example matrices.
matrices to first order, under a certain assumption. We need the following standard
result (see, for example, [21, pp. 183-184]).
Theorem 4.1. Let - be a simple eigenvalue of the matrix A, with right and left
eigenvectors x and y, and let e
a perturbation of A. Then there is an
eigenvalue e
- of e
A such that
We can now prove the following theorem.
Theorem 4.2. Let A 2 C n\Thetan have the Jordan canonical form (2.1), with
j. Suppose that kXk
there is a perturbation e
p-norms, such that
ae( e
Proof. By assumption, - 1 is a simple eigenvalue, so x
are the right eigenvector and the left eigenvector corresponding to - 1 . From Theorem
4.1 we know that any perturbation e
A will have an eigenvalue
Then E has rank 1,
Now for an n \Theta n matrix B and any 1 -
and together with the conditions of the theorem this gives
The proof is completed by taking absolute values in (4.3).
Because of the assumptions on kXk 1 and kX \Gamma1 k1 we do not have the freedom to
choose X to minimize -(X) in Theorem 4.2. We note, however, that if A is diagonalizable
and X has columns all of the same norm, then the condition in Theorem 4.2 on
the rows and columns of X and X \Gamma1 reduces to the requirement that the eigenvalue
of largest modulus be the most ill-conditioned.
Theorem 4.2 enables us to obtain the following corollary of Theorem 3.2.
Corollary 4.3. Suppose that A 2 C n\Thetan is diagonalizable and satisfies the
conditions of Theorem 4.2, and suppose that the O(ffl 2 ) term in (4.2) is negligible. If
is a modest constant depending only on n, then
Proof. By Theorem 4.2, if ae
Rearranging gives
Using Theorem 3.2 we have the required result for c
Suppose we compute the eigenvalues of A by a backward stable algorithm, that
is, one that yields the exact eigenvalues of A
c n a modest constant. (An example of such an algorithm is the QR algorithm [5,
Sec. 7.5]). Then the computed spectral radius b ae satisfies b ae - ae cn ukAk2 (A). In view of
Corollary 4.3 we can formulate a rule of thumb:
The computed powers of A can be expected to converge to zero if the
spectral radius computed via a backward stable eigensolver is less than 1.
This rule of thumb has also been discussed by Trefethen and Trummer [23] and Reichel
and Trefethen [19]. In our experience the rule of thumb is fairly reliable when b
ae is not
too close to 1. For the matrices used in our examples we have
and we observed convergence of the computed powers for C 8 and C 13 + 0:01I and
divergence for the other matrices.
A. Approximating X in the Jordan Form of Cn . In Section 3 we needed an
estimate of -(X) for the Chebyshev differentiation matrix, C, where
the Jordan form. In this appendix we outline our approach for computing an estimate
of -(X).
Recall that
where J is a single Jordan block whose diagonal is zero. Suppose we decompose C
into Schur form via the orthogonal matrix Q (which is real since C's spectrum is real),
that is,
where T is upper triangular with zero diagonal. If we can find an upper triangular
matrix R such that (R). We require
is the jth column
of R.
We choose the arbitrary last column of R to be the last column of the identity
matrix. The following algorithm computes R (here, we use the Matlab colon notation).
It remains to compute the Schur form of C. We do not use the QR algorithm to
compute the Schur form, as for nilpotent matrices it can lead to triangular matrices
with elements of order 1 on the diagonal. We use the following algorithm described by
Golub and Wilkinson in [6, Section 10], which, although computationally expensive,
has good error properties.
Compute the SVD of C
1 .
Compute the SVD C
Upon completion of the algorithm we have lower triangular,
and so we apply our first algorithm to L T to estimate - 2 (X) (note that the Jordan
matrix for A T is a permutation of the one for A [14, Section 3.2.3]).

Acknowledgements

. We thank Des Higham and Nick Trefethen for their comments
on the manuscript.



--R

Eigenvalues of Matrices
Qualitative computing: Elements of a theory for finite precision computation
The Theory of Matrices
The asymptotic behaviour of powers of matrices
Matrix Computations

Bounds for iterates
Stiffness of ODEs
Algorithm 694: A collection of test matrices in MATLAB


Componentwise error analysis for stationary iterative meth- ods

Matrix Analysis

Nonnormality and Jordan condition numbers of matrices
Nineteen dubious ways to compute the exponential of a matrix
Solution of Equations in Euclidean and Banach Spaces
Eigenvalues and pseudo-eigenvalues of Toeplitz matrices
Functions of a matrix and their norms
Matrix Perturbation Theory
Pseudospectra of matrices
An instability phenomenon in spectral methods
Condition numbers and equilibration of matrices
From the Buffon needle problem to the Kreiss matrix theorem
--TR

--CTR
N. D. Assimakis , E. Z. Psarakis , D. G. Lainiotis, Steady state Kalman filter: a new approach, Neural, Parallel & Scientific Computations, v.11 n.4, p.485-490, December
rounding errors;pseudospectrum;jordan canonical form;nonnormal matrices;matrix powers