New Methods for Estimating the Distance to Uncontrollability.
Controllability is a fundamental concept in control theory. Given a linear control system, we present new algorithms for estimating its distance to uncontrollability, i.e., the norm of the normwise smallest perturbation that makes the given system uncontrollable. Many algorithms have been previously proposed to estimate this distance. Our new algorithms are the first that correctly estimate this distance at a cost polynomial in a dimension of the given system. We report results from some numerical experiments that demonstrate the reliability and effectiveness of these new algorithms.
Introduction
One of the most fundamental concepts in control theory is that of controllability. A matrix pair
n\Thetan \Theta C n\Thetam is controllable (see Kailath [20, pages 85-90]) if the state function
in the linear control system
can be directed from any given state to a desired state in finite time by an input
could signal fundamental trouble with the control model or the underlying physical
system itself (Byers [11]).
A large number of algebraic and dynamic characterizations of controllability have been given
(Laub [21], for example). But each and every one of these has difficulties when implemented in
finite precision (Patel, Laub, and Van Dooren [27, page 15]). For instance, it is well known that
(A; B) is controllable if and only if
where C is the set of complex numbers. However, it is not clear how to numerically verify whether a
system is controllable through (1.2). More critically, equation (1.2) does not provide any means to
detect systems that are "nearly" uncontrollable, systems that could be equally troublesome. From
these considerations, it became apparent (see Laub [21] and Paige [26]) that a more meaningful
Department of Mathematics, University of California, Los Angeles, CA 90095-1555. This research was supported
in part by NSF Career Award CCR-9702866 and by Applied Mathematical Sciences Subprogram of the Office of
Energy Research, U.S. Department of Energy under Contract DE-AC03-76SF00098.
measure is the distance to uncontrollability, the norm distance of the pair (A; B) from the set of all
uncontrollable pairs:
It was later shown by Eising [15, 16] that
where oe n (G) denotes the n-th singular value of G 2 C n\Theta(n+m) . Demmel [12] relates ae(A; B) to the
sensitivity of the pole-assignment problem.
Many algorithms have been designed to compute ae(A; B). However, the function to be minimized
in (1.4) is not convex and may have as many as n or more local minima. It is not clear
just how many local minima there are for any given problem (Byers [11]). Methods that search
for a local minimum tend to be efficient but have no guarantee of finding ae(A; B) with any accu-
racy, since ae(A; B) is the global minimum (Boley [4, 6], Boley and Golub [5], Boley and Lu [7],
Byers [11], Elsner and He [17], Miminis [24], and Wicks and DeCarlo [31]); and methods that
search for the global minimum (Byers [11], Gao and Neumann [18], and He [19]) sometimes do have
this guarantee, but require computing time that is inverse proportional to ae 2 (A; B), prohibitively
expensive for nearly uncontrollable systems, the kind of systems for which computing ae(A; B) is
important. While the backward stable algorithms of Beelen and Van Dooren [2, 3, 30] and Demmel
and K-agstr-om [13, 14] are efficient and very useful for detecting uncontrollability, they often fail to
detect near-uncontrollability.
In this paper, we propose new methods to correctly estimate ae(A; B) to within a factor of 2.
They are based on the following bisection method:
Algorithm 1.1 Bisection Method.
while
endwhile
The bisection idea was used to compute the distance of a stable matrix to the unstable matrices
(Byers [10]). It was then used to compute the L1 norm of a transfer matrix (Boyd, Balakrishnan
and Kabamba [9]); a quadratically convergent version of this later method was developed by Boyd
and Balakrishnan [8].
There were past attempts to use Algorithm 1.1 to estimate ae(A; B) as well [11, 18]; but they
have resulted in potentially prohibitively expensive algorithms. The critical difference between our
new approach and earlier attempts lies in how to numerically verify whether ffi - ae(A; B). Our new
approach is based on a novel verifying scheme (see Section 3.2). Paralleling the development of
Boyd and Balakrishnan [8], we have also developed a generally quadratically convergent version of
Algorithm 1.1. With very little modification, our new methods can be used to detect the uncontrollable
modes for any given tolerance. The knowledge of such modes is essential if one wishes to
remove them from the system.
Complexity-wise, these new algorithms differ from previous algorithms in that they are the first
algorithms that correctly estimate the distance at a cost polynomial in the matrix size. In fact,
they require O(n 6 ) floating pointing operations. The main cost of these new algorithms is the
computation of some eigenvalues of certain sparse generalized eigenvalue problems of size O(n 2 ).
In x2 we review methods of Byers and Gao and Neumann to minimize the function in (1.4)
when - is restricted to a straight line on the complex plane. In x3 we present our new methods
to minimize the function in (1.4) over the entire complex plane. In x4 we present some numerical
results. And In x5 we draw conclusions and discuss open questions.
Minimization Methods over a Straight Line
ih
where - is a real variable. To motivate our new methods to minimize the function in (1.4) over
the entire complex plane, in this section we present Algorithm 2.1 below to estimate the global
minimum of g(-) to within a factor of 2 for a given complex number - 0 and a real angle '. This
algorithm is a variation of the bisection schemes of [11, 18], which can actually compute the global
minimum. Let -   be a global minimizer.
Algorithm 2.1 Bisection Method over a Straight Line
while
endwhile
We will also discuss a quadratically convergent version of Algorithm 2.1 in x2.2.
2.1 The Bisection Method over a Straight Line
Let g(-   What is missing in Algorithm 2.1 is a scheme to numerically verify whether
We discuss such a scheme below. Different versions of it were developed
in Byers [11] and Gao and Neumann [18], and were based on earlier work of Byers [10].
We assume that ffi - g(-   ). Since g(-) is a continuous function with
lim
it follows that there exist at least two solutions 1 to the equation By the definition of
singular values, this implies that there exist non-zero vectors
x
y
and z such that
I B
x
y
A   \Gamma
I
x
y
has a double root at -   .
These equations can be rewritten asB B @
I B
A   \Gamma
I \Gammaffi I 0
z
x
To simplify (2.2), we QR-factorize
\Gammaffi I
!/
and define /
z
y
These relations and equation (2.2) imply that R   z must be non-singular. It
follows that z Hence equation (2.2) is reduced to@ A \Gamma
I
\Gammaffi I
A   \Gamma
x
which can be further rewritten as
\Gammaffi I
A
!/
x
e i' I 0
!/
x
Since Q 12 is part of the Q factor in the QR factorization (2.3), it follows that
22
which imply
Hence Q 12 is non-singular and the pencil in (2.4) is regular. It is now easy to show that condition
only if the matrix pencil in (2.4) has a real eigenvalue - .
To verify whether ffi - g(-   ) in Algorithm 2.1, we compute the eigenvalues of the pencil in (2.4).
If this pencil has real eigenvalues, then
guarantees that 2ffi - g(-   ) from the previous bisection step, the value of ffi after Algorithm 2.1 exits
from the while loop must satisfy
We note that equation (2.2) was not reduced in [11], making it more time consuming to verify
whether It was reduced to a regular eigenvalue problem by solving for y in [18], but the
reduction appears to be less numerically reliable than our reduction to (2.4).
2.2 A Quadratically Convergent Variation
Boyd and Balakrishnan [8] note that in the context of computing the L1 -norm of a transfer matrix,
the function to be minimized is in general approximately quadratic near the maximum, and they
used this fact to design a quadratically convergent variation of a bisection method for computing
the L1 -norm.
Their idea applies equally well in minimizing (2.1). If enough and if - 1 and - 2
are two roots of closest to -   , then arguments similar to those of [8] show that (- 1
is in general a much better approximation to -   . We summarize this algorithm below.
Algorithm 2.2 Quadratically Convergent Variation of Algorithm 2.1
while
Choose two real eigenvalues - 1 and - 2 of the pencil (2.4).
ae
/2 .
endwhile
With arguments similar to those of [8], it is easy to show that if - 1 and - 2 are chosen correctly, then
This estimate holds even if g(-) is not approximately quadratic near -   . We caution that strictly
speaking Algorithm 2.2 is not even asymptotically quadratically convergent, since it terminates
as soon as it has found a ffi that satisfies (2.5). Nevertheless, relation (2.6) does indicate rapid
convergence of Algorithm 2.2 when g(-   ) is tiny.
3 Minimization Methods over the Complex Plane
Now we discuss methods to minimize the function in (1.4) over the entire complex plane. Let
One such method is Algorithm 1.1 discussed in x1. As in Algorithm 2.1, we need to develop a
scheme to verify whether ffi ? ae(A; B) in order to complete Algorithm 1.1. To do so, we first prove
a fundamental theorem in x3.1; we then provide such a scheme in x3.2; and finally we develop a
generally quadratically convergent version of Algorithm 1.1 in x3.3.
3.1 A Fundamental Theorem
Our scheme to verify whether ffi ? ae(A; B) is based on Theorem 3.1 below.
Theorem 3.1 Assume that ffi ? ae(A; B). Then there are at least two pairs of real numbers ff and
fi such that
denotes a singular value of G.
Proof: From standard perturbation theory we have
Hence f(ff; fi) goes to infinity if jff does. It is well-known that f(ff; fi) is a continuous function
of ff and fi. Consequently the fact that ffi ? ae(A; B) immediately implies that there exists a pair of
From the definition of singular values, ff and fi satisfy
if and only if they satisfy the algebraic equation
det
It follows from f(ff 1 that this algebraic equation has at least one solution; and it follows
from (3.3) that all its solutions are finite. Consequently, these solutions form a finite number of
closed (continuous) algebraic curves on the ff-fi plane.
Now we claim that the point (ff   ; fi   ) must be in the interior of one of these closed curves. In
fact, if this is not the case, then there exists a continuous curve - 1 (- 2 (-)) on the ff-fi
plane that does not intersect with any of these algebraic curves but "connects" (ff   ; fi   ) and infinity:
In other words,
It follows from the continuity argument that there exists a - 1 such that f (- 1 (- 1
this contradicts the assumption that the curve -) does not intersect with any of the algebraic
curves. Consequently, the point (ff   ; fi   ) must be in the interior of one of these closed curves. Among
all closed curves that have (ff   ; fi   ) in their interior, let G denote the one that covers the smallest
area.
It follows from the same continuity argument that there exist two points
on G with In other words,
This fact has been shown in Byers [11] and Gao and Neumann [18].
For simplicity, we assume that P 1 and P 2 are chosen so that j 1 and j 2 are the smallest positive
numbers.
Since the point (ff   ; fi   ) is in the interior of G and also lies strictly inside the line segment between
that any point that lies strictly inside this line segment is in the interior of
curve G. Combining (3.4) with relation (3.3), we get
Now we shift all the points on G horizontally by the same amount \Gammaj to get a closed curve
is a point on G.g
Since P 2 is a point on G, b
a point on b
G. Assume that
Then relation (3.5) implies that b
is a point that lies strictly inside the line segment between P 1
and P 2 . Hence b
is in the interior of curve G. Let P be the leftmost point on G. Then
is the leftmost point on b
G. we have that b
P 3 is in the exterior of G.
In other words, we have found a point b
that is on b
G and in the interior of G, and another
point b
P 3 on b
G that is in the exterior of G. Since G and b
G are continuous closed curves, we conclude
that these two curves intersect.
intersecting point. It follows that both (ff 4 must be
points on G. Hence ff 4 and fi 4 are a solution to (3.2). Therefore equations (3.2) have at least one
solution.
In the following argument we assume that (ff 4 ; fi 4 ) is the only intersecting point of G and b
G. Let
G 1 denote the set of points on b
G that are either on G or in the interior of G and let G 1 denote the
corresponding set of points on G. If b
G 1 is not a closed curve itself, then b
must be an open curve
with one end point on G and the other in the interior of G. It follows that b
must
have positive arclength. Hence the portion of G without G 1 is a closed curve. But this contradicts
the way G is constructed. This contradiction implies that b
must be closed curves
themselves. Let b
denote the set of points on b
G that are either on G or in the exterior of G and
let G 2 denote the corresponding set of points on G. A similar argument shows that G 2 must be a
closed curve as well.
By construction, b
do not share any common region with positive
area. Hence (ff   ; fi   ) can only be in the interior of one of these closed curves, this implies that G
is not a closed curve that has (ff   ; fi   ) in its interior and covers the smallest area, a contradiction
to the way G was constructed. This contradiction is the result of the assumption that G and b
G
intersect only once. Hence G and b
G must intersect at least twice, so equations (3.2) must have at
least two real solutions. By a continuity argument, equations (3.2) have two, possibly identical, real
solutions, even if
3.2 A New Verifying Scheme
In the following we consider how to numerically verify whether equations (3.2)
have a real solution. By the definition of singular values, equations (3.2) imply that there exist
non-zero vectors
x
y
x
y
, and b
z such that
x
y
x
y
x
x
y
These equations can be rewritten asB @
A   \Gamma ffI \Gammaffi I 0
z
x
z
x
and 0
z
x
z
x
In the QR factorization (2.3), define
z
y
and
b z
These relations and equations (3.6) and (3.7) imply that R   z
non-singular, it follows that z Hence equations (3.6) and (3.7) are reduced to
\Gammaffi I
!/
x
I 0
!/
x
and /
\Gammaffi I
!/
x
I 0
!/
x
As shown in x2.1, Q 12 is always non-singular for Hence the matrix on the right hand sides of
both (3.8) and (3.9) is non-singular. In order for the two pencils defined in (3.8) and (3.9) to share
a common pure imaginary eigenvalue fii, the following matrix equation for X 2 R 2n\Theta2n
\Gammaffi I
I 0
I 0
\Gammaffi I
must have a non-zero solution. Partition
, this matrix equation becomes
12\Omega I 0
vec (X 22 )
vec (X 21 )C C C A ;
12\Omega
)\Omega I 0
I\Omega
\Omega I
12\Omega
12\Omega I \Gammaffi
In these
equations,\Omega is the Kronecker product and vec(G) is a vector formed by stacking the
column vectors of G.
To reduce (3.11) to a standard generalized eigenvalue problem, let
(R
be the RQ factorization of H,
; and define
Then the first equation in (3.11) reduces to setting the second equation in (3.11)
becomes
where
\GammaQ
12\Omega I 0I\Omega Q 12
Equation (3.13) is now a 2n 2 -by-2n 2 generalized eigenvalue problem. Hence we have reduced the
problem of finding a non-zero solution to (3.10) to the generalized eigenvalue problem (3.13).
To summarize, we have shown that in order for (3.2) to have at least one real solution (ff; fi),
both matrix pencils in (3.8) and (3.9) must share a common pure imaginary eigenvalue fii. This
requires that the matrix equation (3.10) must have a non-zero solution, which, in turn, is equivalent
to requiring that the generalized eigenvalue problem (3.13) have a real eigenvalue ff.
In order to verify whether ffi ? ae(A; B) in any bisection step of Algorithm 1.1, we set
in (3.2) and check whether the generalized eigenvalue problem (3.13) has any real eigenvalues ff.
If it does, we then check for each real ff whether the two matrix pencils in (3.8) and (3.9) share a
common pure imaginary eigenvalue fii. If they do for at least one ff, then we have found a pair of
ff and fi such that
On the other hand, if (3.13) does not have a real eigenvalue, or if the matrix pencils in (3.8) and (3.9)
do not share a common pure imaginary eigenvalue for any real eigenvalue of (3.13), then we conclude
by Theorem 3.1 that
On the other hand, Algorithm 1.1 guarantees that 2ffi - ae(A; B) from the previous bisection step.
Thus the value of ffi after Algorithm 1.1 exits from the while loop must satisfy
ae(A; B)
3.3 A Generally Quadratically Convergent Variation
Algorithm 1.1 converges linearly. Following Boyd and Balakrishnan [8] (see x2.2), we develop
a generally quadratically convergent version of Algorithm 1.1 in this section. In the following
development, we assume that f(ff; fi) is analytic in both ff and fi in a small neighborhood of
(ff   ; fi   ) so that f(ff; fi) permits the following expansions
@
@ff
where fl, -, and - are the relevent second-order partial derivatives at (ff   ; fi   ). We further assume
that the matrix \Gamma j
is positive definite. These expansions with a positive definite \Gamma
imply that (ff   ; fi   ) is at least a local minimum.
Now assume that ff and fi are such that f(ff It follows from (3.15) that
@
@ff
Expanding the partial derivative at (ff   ; fi   ) to get
where we have used the fact that @
@ff
be the two solutions to (3.2)
that are near (ff   ; fi   ) (see Theorem 3.1). It follows that ff i and fi i satisfy both (3.17) and (3.16)
2. Now let i 1;i and i 2;i denote the error terms in (3.17) and (3.16) for (ff
Consequently,
It follows from the first equation that
Plugging this into the second equation and simplifying:
where
O
(j
On the other hand, since \Gamma is assumed to be positive definite, equation (3.16) implies that
O(-). Furthermore, the choice of Algorithm 1.1 ensures that
The result can be rewritten as 3
+O(-) and (fi
Combining these equations,
Plugging this relation into (3.18) and combining the equations for
2 and fi new
It follows that
and that
f (ff new ; fi new
We note that this relation is very similar to (2.6). Now we modify Algorithm 1.1 to getThese relations hold as long as flj 2
=4 -. It is likely that under certain conditions, Theorem 3.1 holds for much
larger values of j as well.
Algorithm 3.1 Quadratically Convergent Variation of Algorithm 2.1
while
Choose two real solutions (ff
ae
/2 .
endwhile
In our implementation, we computed f
among adjacent pairs of real solutions
and chose the pair with smallest f value. We note that in both Algorithms 1.1 and 3.1, we
can compute a better initial guess ffi by using Algorithm 2.2 with some values of - 0 and ', such as
Algorithm 3.1 was derived under the assumptions at the beginning of x3.3, which need not
hold for all linear control systems of the form (1.1). Hence estimate (3.19) may not hold for
some linear control systems. However, it is clear that Algorithm 3.1 converges at least linearly to
arrive at an estimate ffi that satisfies (3.14). Similar to Algorithm 2.2, Algorithm 3.1 is not strictly
speaking quadratically convergent since it terminates as soon as it has found a ffi that satisfies (3.14).
Nevertheless, Algorithm 3.1 does converge much more rapidly than Algorithm 1.1 when ae(A; B) is
tiny. We discuss this point further in x4.
3.4 Further Considerations
Sometimes it may be more important to find the uncontrollable modes of (1.1) for a given tolerance
". In this case, we solve equations (3.2) with ". If there are no solutions to (3.2), then the
system (1.1) is controllable; otherwise, each solution to (3.2) corresponds to an uncontrollable mode.
Conversely, it is easy to see from the proof of Theorem 3.1 that any uncontrollable mode will result
in at least two solutions to (3.2). Hence the set of all solutions to (3.2) provide approximations
to the uncontrollable modes of (3.2). The formulas for ff new and fi new provide more accurate
approximations to these modes.
If ae(A; B) is very small, then small during the execution of Algorithms
1.1 and 3.1. In fact, for small enough j, the two different points ff + fii and ff
look identical. Hence the solutions to (3.2) are potentially ill-conditioned. See x4 for more details.
Like many other algorithms in engineering computations, such as those for semi-definite programming
[1, 25, 29], both Algorithms 1.1 and 3.1 are expensive for large problems, since both the
reduction to and the solution of the pencil (3.13) require O(n 6 ) floating point operations. However,
the eigenvalue problem (3.11) is highly sparse as a 4n 2 \Theta 4n 2 problem. It is likely that sparse matrix
computation technologies, such as the implicitly restarted Arnoldi iteration [22, 23, 28], can be used
to compute the real eigenvalues of (3.13) quickly. The effectiveness of this approach is currently
under thorough investigation.
4 Numerical Experiments
We have done some elementary numerical experiments with Algorithms 1.1 and 3.1. In this section
we report some of the results obtained from these experiments. The experiments were done in
matlab in double precision.
Matrices in Examples 2 through 5 were taken from Gao and Neumann [18]. These are systems
with small ae(A; B). Global optimization methods (Byers [11], Gao and Neumann [18], and He [19])
could require prohibitively expensive computation time to correctly estimate ae(A; B) in these cases.
On the other hand, both Algorithms 1.1 and 3.1 worked well on them, with Algorithm 3.1 converging
much faster than Algorithm 1.1 as expected.
Example 1. In this example we took A 2 R 5\Theta5 and B 2 R 5 to be random matrices. This is a
matrix pair with fairly large ae(A; B). Both Algorithms 1.1 and 3.1 took 2 iterations to terminate.
This example illustrates that for linear systems (1.1) that are far away from the set of uncontrollable
systems, both Algorithms 1.1 and 3.1 take very few iterations.
Example 2. In this example we took
This pair is uncontrollable since the smallest singular value of [A \Gamma (1 \Sigma 2i) I B] is zero. Algorithms
1.1 and 3.1 took 42 and 5 iterations, respectively, to find ae(A;
Example 3. In this example we took
\Gamma0:32616458 \Gamma0:09430266 0:05207847 \Gamma:08481401 0:05829280
0:01158922 \Gamma:39787419 \Gamma:14901699 \Gamma:01394125 \Gamma:10626942
0:05623810 \Gamma:03153954 \Gamma:50160557 \Gamma:05748511 \Gamma:00552321
iterations to return returned six distinct
solutions to (3.2): (ff;
and
On the other hand, Algorithm 3.1 took 4 iterations to return
returned one distinct solution (ff;
Example 4. In this example we took
\Gamma:22907968 0:08886286 \Gamma:18085425 \Gamma:03469234 \Gamma:32819211
\Gamma:02507663 :30736050 \Gamma:24819024 :21852948 \Gamma:06260819
iterations
to return returned four distinct
solutions
On the other hand, Algorithm 3.1 took 3 iterations to return
returned two distinct solutions (ff;
Example 5. In this example we took
\Gamma:27422658 \Gamma:21968089 \Gamma:21065336 \Gamma:22134064 0:19235875
\Gamma:07210867 :18848014 \Gamma:29068998 :28936270 0:10007703
\Gamma:03547166 :17931676 :14590007 :00556579 :38838791
\Gamma:07780546 \Gamma:29477373 :01366200 :32749991 \Gamma:0131683C C C C C A
iterations
to return returned 14 distinct solutions
to (3.2). On the other hand, Algorithm 3.1 took 2 iterations to return
for returned one distinct solution (ff;
Example 6. In this example we took the matrix pair (A; B) from Example 2, and set
where Q is a random orthogonal matrix. This new matrix pair is still uncontrollable. But Algorithm
1.1 took 28 iterations to return iterations
to return This example illustrates that both algorithms can have numerical
difficulties in correctly estimating ae(A; B) if it is very tiny.
5 Conclusions and Extensions
In this paper, we have presented the first algorithms that require a cost polynomial in the matrix
size to correctly estimate the controllability distance ae(A; B) for a given linear control system. And
we have demonstrated their effectiveness and reliability through some numerical experiments.
The biggest open question is how to further reduce the cost. At the core of these algorithms
is the computation of all real eigenvalues of a sparse 4n 2 \Theta 4n 2 eigenvalue problem. Currently, we
find these eigenvalues by treating the eigenvalue problem as a dense one, resulting in algorithms
that are too expensive for large problems. In the future, we plan to exploit the possibility of finding
these real eigenvalues via sparse matrix computation technologies, such as the implicitly restarted
Arnoldi iteration [22, 23, 28], to significantly reduce the computation cost;
Another open question is to better understand the effects of finite precision arithmetic on the
estimated distance ae(A; B). As we observed in x4, if ae(A; B) is very tiny, then the distance estimated
by the new algorithms in finite precision could be much larger than the exact distance.
Finally, the perturbation [\DeltaA; \DeltaB] in (1.3) can be complex even if both A and B are real. It is
known (Byers [11]) that the norm-wise smallest real perturbation can be much larger than ae(A; B).
Whether our new algorithms shed new light on the computation of the norm-wise smallest real
perturbation remains to be seen.



--R


An improved algorithm for the computation of Kronecker's canonical form of a singular pencil.
A class of staircase algorithms for generalized state space systems.
Computing the controllability/observability decomposition of a linear time-invariant dynamic system: a numerical approach
The Lanczos-Arnoldi algorithm and controllability
Computing rank-deficiency of rectangular matrix pencils
Measuring how far a controllable system is from uncontrollable one.
A regularity result for the singular values of a transfer matrix and a quadratically convergent algorithm for computing its L1
A bisection method for computing the H1 norm of a transfer matrix and related problems.
A bisection method for measuring the distance of a stable matrix to the unstable matrices.
Detecting nearly uncontrollable pairs.
On condition numbers and the distance to the nearest ill-posed problem


The distance between a system and the set of uncontrollable systems.
Between controllable and uncontrollable.
An algorithm for computing the distance to uncontrollability.
A global minimum search algorithm for estimating the distance to uncontrollability.
Estimating the distance to uncontrollability: A fast method and a slow one.
Linear Systems.
Survey of computational methods in control theory.
Analysis and Implementation of an Implicitly Restarted Arnoldi Iteration.
Deflation techniques for an implicitly restarted Arnoldi iteration.
Numerical algorithms for controllability and eigenvalue location

Properties of numerical algorithms related to computing controllability.
Numerical Linear Algebra Techniques For Systems and Control.
Implicit application of polynomial filters in a k-step Arnoldi method

The computation of Kronecker's canonical form of a singular pencil.
Computing the distance to an uncontrollable system.
--TR
QR factorization;controllability;complexity;numerical stability