SOLVING SINGULAR GENERALIZED EIGENVALUE PROBLEMS
BY A RANK-COMPLETING PERTURBATION
MICHIEL E. HOCHSTENBACH∗, CHRISTIAN MEHL†,AND BOR PLESTENJAK‡
Abstract. Generalized eigenvalue problems involving a singular pencil are very challenging to
solve, both with respect to accuracy and efficiency. The existing package Guptri is very elegant but
may sometimes be time-demanding, even for small and medium-sized matrices. We propose a simple
method to compute the eigenvalues of singular pencils, based on one perturbation of the original
problem of a certain specific rank. For many problems, the method is both fast and robust. This
approach may be seen as a welcome alternative to staircase methods.
Key words. Singular pencil, singular generalized eigenvalue problem, rank-completing pertur-
bation, Guptri, model updating, double eigenvalues, two-parameter eigenvalue problem, differential
algebraic equations, quadratic two-parameter eigenvalue problem.
AMS subject classifications. 65F15, 15A18, 15A22, 15A21, 47A55, 65F22
1. Introduction. We study the computation of eigenvalues of small to medium-
sized matrix pencils A−λB, where Aand Bare (real or complex) n×mmatrices and
where, in addition, the matrix pencil A−λB is singular, which means that m6=nor
if m=nthen
det(A−λB)≡0.
In these cases, the common definition of eigenvalues as roots of det(A−λB) obviously
does not make sense. Therefore, finite eigenvalues of a singular matrix pencil A−λB
are typically defined as values λ0∈Csatisfying rank(A−λ0B)<nrank(A, B), where
nrank(A, B) := max
ζ∈Crank(A−ζB)
denotes the normal rank of the pencil A−λB; see [10]. Similarly, we say that ∞is an
eigenvalue of the singular pencil A−λB if rank(B)<nrank(A, B). In the following
we will restrict ourselves to the case m=nas the case m6=ncan easily be reduced
to the square case by adding an appropriate number of zero rows or columns.
The singular generalized eigenvalue problem (singular GEP) is well known to be
ill-conditioned as arbitrarily small perturbation may cause drastic changes in the eigen-
values. A classical example is given by the pencils
A−λB =1 0
0 0 −λ1 0
0 0 and e
A−λe
B=1ε1
ε20−λ1ε3
ε40,
∗Version May 22, 2018. Department of Mathematics and Computer Science, TU Eindhoven, PO
Box 513, 5600 MB, The Netherlands, www.win.tue.nl/∼hochsten. This author has been supported
by an NWO Vidi research grant.
†Institut f¨ur Mathematik, Technische Universit¨at Berlin, Straße des 17. Juni 136, 10623 Berlin,
grant.
‡IMFM and Faculty of Mathematics and Physics, University of Ljubljana, Jadranska 19, SI-1000
Slovenian Research Agency (grant P1-0294).
1
2HOCHSTENBACH, MEHL, AND PLESTENJAK
where ε1, . . . , ε4∈C\ {0}; see [22]. While A−λB is singular and has only the
eigenvalue 1, the perturbed pencil e
A−λe
Bis regular and has the eigenvalues ε1
ε3and ε2
ε4
that can be anywhere in the complex plane even for tiny absolute values of ε1, . . . , ε4.
On the other hand, it was observed in [43] that a situation as above is exceptional
and that generically small perturbations of a singular square pencil make the pencil
regular and some of the eigenvalues of the perturbed pencil are very close to the
original eigenvalues of the singular pencil. The following example illustrates this. The
Matlab commands
A = diag([1 2 3 0 0 0]);
B = diag([2 3 4 0 0 0]);
eig(U’*A*V, U’*B*V)
where Uand Vare random 6 ×6 orthogonal matrices, give the eigenvalues
0.5000 0.6667 0.7500 0.1595 0.6756 0.6543
We see that the three (finite) eigenvalues of the regular part are correct. Following the
terminology of [40], the other three values are “fake eigenvalues” and correspond to
the singular part of the pencil. (Explicit error analysis for the regular eigenvalues of
singular pencils has been undertaken in [9, 11].) Despite this observation, Van Dooren
suggests in [40] to solve the singular generalized eigenvalue problem by first extracting
the regular part and then use the QZ algorithm on that part. Wilkinson strongly
supports that recommendation in [43].
A robust software package which follows Van Dooren’s recommendation is Gup-
tri [13, 17]. If the pencil is singular, then first a “staircase” algorithm is applied to
deflate the singular part of the pencil, and then the QZ algorithm is used to compute
the eigenvalues of the remaining regular part. While the results of Guptri are usu-
ally excellent, this method may turn out to be time-consuming; for instance, applying
Guptri on a singular 300×300 pencil on our machine took over 20 seconds, while Mat-
lab’s eig on a random pencil of the same size took less than a second.1Another issue
is the fact that the rank decisions, needed in staircase type methods, are involved. If
the pencil has a minimal index of size η, then at least η+ 1 such decisions have to be
taken. Typically, these decisions tend to become more and more critical during the
run of the staircase algorithm. See, e.g., [14] or [33, Ex. 18], where a variant of the
staircase algorithm for the singular two-parameter eigenvalue problem, introduced in
[31], fails in double precision but gets the right result in higher precision.
Another way of extracting the regular part that needs fewer rank decisions was
suggested in [30]. One may interpret the singular pencil as a constant coefficient
differential-algebraic equation and perform a regularization procedure with the help
of an derivative array as described in [3]. In this way, the regular part of the pencil
can be extracted by performing only three nullspace computations. However, the
derivative array approach leads to an inflation of the system by a factor of at least
η+ 1, where ηis the largest minimal index of the given pencil, and may therefore
result in high computational costs.
In this paper, we propose a new method to compute the regular eigenvalues of a
singular pencil. The method is based on considering perturbations of rank
k=n−nrank(A, B)
1We note that this experiment has been performed some years ago. A current practical issue is
that there is no publicly available 64-bit Guptri code.
SOLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION 3
that we will call rank-completing perturbations as the rank is exactly large enough to
generically turn the pencil into a pencil of full normal rank. As we will show, the
canonical form of the original regular part of the given pencil stays invariant under
generic rank-completing perturbations.
The idea of computing eigenvalues of singular pencils with rank-completing pertur-
bations is not completely new, and the following specific type has been used in system
theory as early as in the 70s (without the use of the terminology “rank-completing
perturbation”). If a system of the form
˙x=Ax +Bu
y=Cx +Du
is given, where A∈Rn,n,B∈Rn,m,C∈Rr,n, and D∈Rr,m are the system matrices,
xstands for the state of the system, uis the input, and yis the output, then the
eigenvalues of the system pencil
S(λ) = λI −A B
−C D
are of particular interest in control theory; see [15] and the references therein. (If the
system is minimal, then these eigenvalues are also referred to as transmission zeros of
the system.) Clearly, if m6=r, then the pencil S(λ) is rectangular and thus singular.
For that case and under the additional assumptions r < m and nrank S(λ) = n+r,
the following algorithm which is based on ideas of [6] has been proposed in [24] for
the computation of the transmission zeros:
1: Select two random matrices [E1F1],[E2F2]∈Rm−r,n+mso that
S1(λ) :=
λI −A B
−C D
E1F1
and S2(λ) :=
λI −A B
−C D
E2F2
are regular.
2: Compute the eigenvalues Eiof Si(λ) for i= 1,2.
3: Compute the intersection E=E1∩ E2.
Since for each regular eigenvalue λ0of S(λ) we obviously have nrank Si(λ0)< n+m
for i= 1,2, it immediately follows that the eigenvalues of S(λ) are contained both in
the spectrum of S1(λ) and in the spectrum of S2(λ). Both extended matrix pencils
will also give rise to two sets of fake eigenvalues. Since generically these sets will be
disjoint if the applied perturbations are generated randomly, it follows that the set E
will generically coincide with the set of regular eigenvalues of S(λ).
As pointed out in [15], this method may encounter difficulties in distinguishing
the finite zeros from the infinite ones, in particular if the latter occur with a high
multiplicity. Another problem may occur in identifying which values belong to the
intersection Eand which do not. Although the original eigenvalues of the pencil theo-
retically coincide with a subset of each of the sets of eigenvalues of the two perturbed
pencils, they may still differ slightly in practice due to finite precision arithmetic.
Therefore, a tolerance has to be prescribed that decides when two values are consid-
ered to be equal. If this tolerance is chosen to be too small, then some of the regular
eigenvalues may be missed. If, on the other hand the tolerance is set too high, then
two fake eigenvalues of S1(λ) and S2(λ) that happen to be close are likely to be falsely
identified as a regular eigenvalue.
4HOCHSTENBACH, MEHL, AND PLESTENJAK
In this paper, we show that the regular eigenvalues of a singular pencil can be
efficiently computed with the help of just one rank-completing perturbation of the
form
A−λB +τ U(DA−λDB)V∗,
where U,Vare n×kmatrices with orthonormal columns, DA, DBare diagonal k×k
matrices, and τa nonzero real scalar. The orthonormality of the columns is not strictly
necessary, but convenient, for instance since the norm of the perturbation can easily be
controlled in this case with the help of the parameter τ. The problem of identifying the
subset of regular eigenvalues from the computed eigenvalues of the perturbed pencil
is then taken care of by the key observation that the left and right eigenvectors that
correspond to the regular eigenvalues satisfy orthogonality relations with respect to the
matrices Uand V. Thus, instead of comparing two spectra of two different pencils the
regular eigenvalues can be separated from the fake eigenvalues by using information
from the corresponding left and right eigenvectors. We note that perturbations of
singular matrix pencils have already been considered in [4, 9, 27, 38, 39], but it seems
that a detailed investigation of rank-completing perturbations is new, except for [28],
where the case of singular Hermitian pencils of normal rank n−1 was considered.
The rest of this paper is organized as follows. In Section 3 we review some motivat-
ing applications where one is interested in computing eigenvalues of a singular matrix
pencil. The main theoretical results are presented in Section 4, while the method,
which is based on these results, is introduced in Section 5, followed by some numerical
experiments in Section 6. In Section 7 we discuss singular two-parameter eigenvalue
problems and present a new numerical method for such problems. We summarize
some conclusions in Section 8.
2. Preliminaries. Throughout the paper, we will interpret matrix pencils both
as pairs of matrices (A, B)∈Cn,m ×Cn,m or as n×mmatrix polynomials A−λB of
degree at most one and we will switch between these notation whenever it is useful.
An important tool in the theory of singular pencils is the Kronecker canonical form
(KCF) of a pencil A−λB; see, e.g., [16].
Theorem 2.1 (Kronecker canonical form). Let A−λB be a complex n×mmatrix
pencil. Then there exist nonsingular matrices P∈Cn,n and Q∈Cm,m such that
(2.1) P(A−λB)Q=R(λ) 0
0S(λ),R(λ) = J−λIr0
0Is−λN
with J, N in Jordan normal form and in addition Nbeing nilpotent, and
S(λ) = diagLm1(λ), . . . , Lmk(λ), Ln1(λ)>, . . . , Ln`(λ)>,
where Lj(λ) = [0 Ij]−λ[Ij0] is of size j×(j+ 1) and mi≥0for i= 1, . . . , k and
ni≥0for i= 1, . . . , `.
The pencils R(λ) and S(λ) in Theorem 2.1 are called the regular and the singular
part of A−λB, respectively. The eigenvalues of Jare exactly the finite eigenvalues of
A−λB, while the eigenvalue 0 of Ncorresponds to the infinite eigenvalue of A−λB.
The parameters m1, . . . , mkare called the right minimal indices and the parameters
n1, . . . , n`are called the left minimal indices of A−λB. One can easily check that
the normal rank nrank(A, B) of A−λB is then given by min(n−`, m −k). In the
SOLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION 5
remainder of the paper, we will consider square pencils A−λB, i.e., we have n=m.
Note that this implies k=`, i.e., we have must have the same number of right and left
minimal indices. However, the particular values of the left and right minimal indices
may be distinct.
In contrast to the eigenvalues of singular pencils, the corresponding eigenvectors
and deflating subspaces are not well defined. To understand why, we consider the
following example borrowed from [29] (and slightly adapted). The pencil
(2.2) A−λB =
1−λ0 0
0−λ1
0 0 0
obviously has the regular part R(λ) = [1−λ] and thus the pencil A−λB has the single
eigenvalue λ0= 1 with algebraic multiplicity one. Nevertheless, any vector of the form
x(α, β) := [α β β]>with α, β ∈Csatisfies λ0Bx =Ax and thus could be interpreted
as an eigenvector of the pencil. One may argue that in this example the choice α6= 0
and β= 0 seems to be canonical and gives rise to a unique one-dimensional deflating
subspace “corresponding” to the regular part of the pencil. But on the other hand it
follows from the equality
1/α 0 0
−β/α 1 0
0 0 1
1−λ0 0
0−λ1
0 0 0
α0 0
β1 0
β0 1
=
1−λ0 0
0−λ1
0 0 0
that for any choice of α, β with α6= 0 also the vector x(α, β) can be used to extract the
regular part of the pencil (and thus could also be considered as “corresponding” to the
regular part). For this reason, we restrict ourselves to the computation of eigenvalues
of singular pencils, but do not consider corresponding eigenvectors.
Instead of eigenvectors and deflating subspaces, the concept of reducing subspaces
that was introduced in [41] is more adequate in the case of singular pencils. We say
that a subspace Mis a reducing subspace for the pencil A−λB if dim(AM+BM) =
dim(M)−k, where kis the number of right singular blocks. In the example above, the
reducing subspace associated with the eigenvalue λ0= 1 is exactly given by all vectors
x(α, β) with α, β ∈C. The minimal reducing subspace MRS(A, B) is the intersection
of all reducing subspaces. It is unique and can be numerically computed in a stable
way from the generalized upper triangular form (Guptri), see, e.g., [13]. Guptri exists
for every pencil A−λB and has the form
P∗(A−λB)Q=
Ar−λBr× ×
0Areg −λBreg ×
0 0 Al−λBl
,
where matrices Pand Qare unitary, Ar−λBrhas only right singular blocks in its
KCF, Al−λBlhas only left singular blocks in its KCF, and Areg −λBreg has only
regular blocks in its KCF.
3. Motivation and applications. Before proposing our new method, we first
review some motivating applications besides the computation of transmission zeros
mentioned in the introduction where one wants to compute the eigenvalues of a singular
pencil.
6HOCHSTENBACH, MEHL, AND PLESTENJAK
3.1. Differential algebraic equations and descriptor systems. Linear dif-
ferential algebraic equations (DAEs) with constant coefficients have the general form
E˙x=Ax +f(t), x(t0) = x0,
where E, A ∈Rk,n,t0∈R,x0∈Rn, and f: [t0,∞)→Rkis a given inhomogeneity;
see [23]. Linear time-invariant descriptor systems consist of a DAE combined with a
system input and output and take form
E˙x=Ax +Bu, x(t0) = x0,
y=Cx +Du.
where, in addition, B∈Rk,m,C∈Rp,n,D∈Rp,m. Here, x: [t0,∞)→Rnstands for
the state of the descriptor system, u: [t0,∞)→Rmis the input and y: [t0,∞)→Rp
the output. As is highlighted in [3], the problem may be well-posed even if the under-
lying pencil A−λE is singular. Indeed, even in the singular case the corresponding
DAE may have a (even unique) solution for particular inhomogeneities fand for
special initial conditions (t0, x0).
3.2. Double eigenvalue problem. Given two n×nmatrices Aand B, we are
interested in all values λsuch that A+λB has a double eigenvalue. For the generic
case of a double eigenvalue, we look for independent vectors xand ysuch that
(A+λB −µI)x= 0,
(A+λB −µI)2y= 0.
This application is discussed in [33], together with a staircase algorithm to solve it;
see also [21]. In the generic case, the problem has n(n−1) solutions.
In [33], the problem is solved by linearizing the second equation first and then
solving the obtained singular two-parameter eigenvalue problem. The values λthat
we are looking for, are eigenvalues of the singular pencil ∆1−λ∆0of size 3n2×3n2,
where
∆1=A⊗R−I⊗P, ∆0=B⊗R−I⊗Q,
for
P=
A2AB +BA −2A
0I0
0 0 I
, Q =
0B2−B
−I0 0
0 0 0
, R =
0−B I
0 0 0
−I0 0
.
For more details, see [33] as well as [19] for other possible linearizations.
3.3. Singular two-parameter eigenvalue problems. In the case of a singular
two-parameter eigenvalue problem (2EP) one has a pair of singular pencils ∆1−λ∆0
and ∆2−µ∆0(see also (7.2) and (7.3)), and the goal is to find finite regular eigenvalues
(λ0, µ0), where (see, e.g., [31] for more details):
•λ0is a finite regular eigenvalue of ∆1−λ∆0,
•µ0is a finite regular eigenvalue of ∆2−µ∆0,
•there exists a common regular eigenvector z6= 0 such that (∆1−λ0∆0)z= 0,
(∆2−µ0∆0)z= 0, and z6∈ MRS(∆i,∆0) for i= 1,2.
This problem requires more than just solving one singular GEP. We discuss it in more
details in Section 7 and present a new numerical method for its solution.
SOLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION 7
4. Rank-completing perturbations of singular pencils. Let A−λB be a
singular pencil, where A, B ∈Cn,n and nrank(A, B) = n−k. In this section we inves-
tigate the effect of rank-completing perturbations, i.e., rank-kgeneric perturbations of
the form
(4.1) e
A−λe
B:= A−λB +τ(UDAV∗−λ UDBV∗),
where DA, DB∈Ck,k are diagonal matrices such that DA−λDBis a regular pen-
cil, U, V ∈Cn,k have full column rank, and τ∈Ris nonzero. We investigate the
dependence of eigenvalues and eigenvectors of the perturbed pencils on τ.
Above and in the following, the term generic is understood in the following sense:
a set A ⊆ Cmis called algebraic, if it is the set of common zeros of finitely many
complex polynomials in mvariables, and Ais called proper if A 6=Cm. A set Ω ⊆Cm
is called generic if its complement is contained in a proper algebraic set. In this sense
we say that a property Pholds generically with respect to the entries of Uand V,
if there exists a generic set Ω ⊆(Cn,k)2(whereby we interpret (Cn,k)2as C2nk) such
that Pholds for all pencils of the form (4.1) with (U, V )∈Ω.
Remark 4.1. Perturbations of the form (4.1) are not the most general perturba-
tions of rank k. Indeed, completing Uand Vto nonsingular matrices P= [Ue
U] and
Q= [Ve
V], we obtain that
P−1(UDAV∗−λUDBV∗)(Q∗)−1=DA−λDB0
0 0 ,
which means that the perturbation pencil UDAV∗−λUDBV∗has the regular part
DA−λDBof size k×k, while, generically, a matrix pencil of rank k < n would have
no regular part [8]. Thus, more generally one could consider perturbations of the form
(4.2) (A+U1V∗
1+V2U∗
2, B +U1W∗
1+W2U∗
2),
where `∈ {0,1, . . . , k}and where U1, V1, W1∈Cn,` and U2, V2, W2∈Cn,k−`have full
column rank. However, we will restrict ourselves to perturbations of the form (4.1),
because of their favorable properties.
Generically, a rank completing perturbation of A−λB will result in a pencil of full
normal rank, i.e., the perturbed pencil is regular. We will show in the following, that if
rank-completing perturbations of the form e
A−λe
B:= A−λB+τ(UDAV∗−λUDBV∗)
as in (4.1) are considered, then generically the KCF of e
A−λe
Bis for all τ∈R\ {0}
given by
Rreg(λ) 0 0
0Rpre(λ) 0
0 0 Rran(λ)
,
where Rreg(λ) is the regular part of the original pencil A−λB,Rpre(λ) = DA−λDB
and Rran(λ) only has simple eigenvalues that are different from the eigenvalues of
Rreg(λ) and Rpre(λ). Thus, the eigenvalues of e
A−λe
Bare exactly the pregular eigen-
values of A−λB (counted with multiplicities) and n−pnewly generated eigenvalues
which consist of k“prescribed” eigenvalues which are the eigenvalues of the perturba-
tion pencil U(DA−λDB)V∗, and n−p−k“random” eigenvalues that are gathered
in the part Rran(λ).
8HOCHSTENBACH, MEHL, AND PLESTENJAK
We start by showing that under a rank-completing perturbation, the regular part
of A−λB will stay invariant in the above sense. Although our main focus are square
pencils, we will state some results in more generality covering also the case of rect-
angular pencils. The following proposition is a generalization of Theorem 4.2 in [28]
(which dealt with Hermitian pencils) to a block case without any specific structure in
matrices Aand B.
Proposition 4.2. Let A−λB be an n×msingular matrix pencil having at least
kleft minimal indices. Furthermore, let U∈Cn,k. Then generically (with respect to
the entries of U) there exist nonsingular matrices P, Q such that
P(A−λB)Q=R(λ) 0
0S(λ)and PU =0
e
U,
where R(λ)and S(λ)are the regular and singular parts of A−λB, respectively, and
PU is partitioned conformably with P(A−λB)Q.
Proof. Without loss of generality we may assume that A−λB is already in the
KCF
A−λB =R(λ) 0
0S(λ),R(λ) = J−λI 0
0I−λN ,
where S(λ) is the singular part and R(λ) is the regular part of A−λB with J, N
being in Jordan normal canonical form and Nbeing in addition nilpotent. Our main
strategy is to use the part of Uthat corresponds to karbitrarily chosen left minimal
indices to introduce zeros in the components of Uthat correspond to the regular part
of A−λB. Here, we can treat the components of Ucorresponding to the “finite
eigenvalue part” J−λI and the “infinite eigenvalue part” I−λN separately and
we will give the proof only for the first case as the proof for the “infinite eigenvalue
part” is completely analogously. Since we will only transform parts of the singular
part S(λ) that correspond to karbitrarily chosen left minimal indices and leave all
other parts of S(λ) unchanged, it is sufficient to assume that S(λ) consists of only k
singular blocks corresponding to kleft minimal indices and that R(λ) does not have a
part corresponding to infinite eigenvalues. This assumption will simplify the notation
considerably.
Thus, we assume that Uand A−λB have the forms
u=
k
rU0
n1+1 U1
.
.
..
.
.
nk+1 Uk
and A−λB =
r n1. . . nk
rJ−λI
n1+1 Ln1(λ)>
.
.
....
nk+1 Lnk(λ)>
.
In the following, we will use the notation G(d)
`= [0 I`]>and G(u)
`= [I`0]>such
that we can write L`(λ)>=G(d)
`−λG(u)
`for a singular block corresponding to a left
minimal index `. For the transformation matrices P, Q we make the ansatz
P=
r n1+1 . . . nk+1
rI P1. . . Pk
n1+1 I
.
.
....
nk+1 I
and Q=
r n1. . . nk
rI−P1G(u)
n1. . . −PkG(u)
nk
n1I
.
.
....
nkI
SOLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION 9
Since Pand Qare designed in such a way that Bremains unchanged, we obtain
P(A−λB)Q=
J−λI P1G(d)
n1−JP1G(u)
n1. . . PkG(d)
nk−JPkG(u)
nk
Ln1(λ)>
...
Lnk(λ)>
.
It remains to find solutions Pito the equations PiG(d)
ni−JPiG(u)
ni= 0 for i= 1, . . . , k
to obtain P(A−λB)Q=A−λB. Setting Pi= [pi,0pi,1. . . pi,ni] for i= 1, . . . , k,
the equations to be solved take the form
[pi,1. . . pi,ni] = PiG(d)
ni=JPiG(u)
ni= [Jpi,0. . . Jpi,ni−1], i = 1, . . . , k.
Thus, we may choose Pi= [pi,0Jpi,0. . . Jnipi,0], where pi,0∈Cris arbitrary. We
will now use the freedom in the choice of pi,0to guarantee that PU has the desired
form. To this end, let ui,0, . . . , ui,ni∈Ckbe the rows of Ui,i= 1, . . . , k. Then the first
block component of P U is given by U0+P1U1+· · · +PkUkand to set this component
to zero, we have to solve the equation
(4.3) −U0=
k
X
i=1
PiUi=
k
X
i=1
ni
X
j=0
Jjpj,0u>
i,j
for p1,0, . . . , pk,0. Using the vec-operation that “vectorizes” a matrix by stacking the
columns on top of another and recalling the well-known identity vec(XY Z) = (Z>⊗
X) vec(Y) for matrices X, Y, Z, where ⊗denotes the Kronecker product, we obtain
that
(4.4) −vec(U0) =
k
X
i=1
ni
X
j=0 ui,j ⊗Jjpi,0=M·p1,0
.
.
.
pk,0,
where
M=hn1
X
j=0
u1,j ⊗Jj. . .
nk
X
j=0
uk,j ⊗Jji∈Crk,rk.
Clearly, the determinant of Mis a polynomial in the nk entries of U(in fact, it
only depends on the entries of U1, . . . , Uk) which is nonzero for the particular choice
u1,0=e1, . . . , uk,0=ekand ui,j = 0 for j > 0, where e1, . . . , ekdenote the standard
basis vectors of Ck. (Indeed, in this case Mis just the identity of size rk ×rk.)
Thus, generically (with respect to the entries of U), the matrix Mis invertible, so
the equation (4.4) and thus also (4.3) can be uniquely solved for p1,0, . . . , pk,0which
finishes the proof.
The next result shows that the canonical form of the regular part of the original
pencil stays invariant under a generic rank-completing perturbation of the form (4.1).
Concerning the eigenvalues of the perturbed pencil that are also eigenvalues of the
original singular pencil, the result also states that the corresponding left and right
eigenvectors satisfy a particular orthogonality relation.
Theorem 4.3. Let A−λB be an n×nsingular pencil of normal rank n−k, let
U, V ∈Cn,k have full column rank and let DA, DB∈Ck,k be such that DA−λDBis
10 HOCHSTENBACH, MEHL, AND PLESTENJAK
regular and all eigenvalues of DA−λDBare distinct from the eigenvalues of A−λB.
Then, generically with respect to the entries of Uand V∗, the following statements
hold for the pencil (4.1):
1. For each τ6= 0, there exists nonsingular matrices e
Pand e
Qsuch that
(4.5) e
Pe
A−λe
Be
Q=R(λ) 0
0Rnew(λ),
where R(λ)is the regular part of the original pencil A−λB, and Rnew(λ)is
regular and all its eigenvalues are distinct from the eigenvalues of R(λ).
2. If λ0is a finite regular eigenvalue of A−λB, i.e., rank(A−λ0B)< n −k,
then λ0is an eigenvalue of (4.1) for each τ6= 0. Furthermore, the right null
space
Nr(λ0) := ker A+τ UDAV∗−λ(B+τ UDBV∗)
and the left null space
Nl(λ) := ker A+τ UDAV∗−λ(B+τ UDBV∗)∗
are both constant in τ6= 0. In addition:
(a) Nr(λ0)⊥span(V), i.e., if xis a right eigenvector of e
A−λe
Bassociated
with λ0, then V∗x= 0.
(b) Nl(λ0)⊥span(U), i.e., if yis a left eigenvector of e
A−λe
Bassociated with
λ0, then U∗y= 0.
3. If ∞is a regular eigenvalue of A−λB, i.e., rank(B)< n −k, then ∞is an
eigenvalue of (4.1) for each τ6= 0. The right and left null spaces
Nr(∞) := ker(B+τ UDBV∗)and Nl(∞) := ker (B+τ UDBV∗)∗
are both constant in τ6= 0. In addition:
(a) Nr(∞)⊥span(V), i.e., if xis a right eigenvector of e
A−λe
Bassociated
with ∞, then V∗x= 0.
(b) Nl(∞)⊥span(U), i.e., if yis a left eigenvector of e
A−λe
Bassociated with
∞, then U∗y= 0.
Proof. First, we will assume that A−λB does not have one of the eigenvalues 0
or ∞and we will show 1) and 2) for this particular case.
Applying Proposition 4.2, there exist nonsingular matrices P, Q such that
(4.6) P(A−λB)Q=R(λ) 0
0S(λ), PU =0
U2, Q∗V=V1
V2,
where R(λ) and S(λ) are the regular and singular parts of A−λB, respectively, both
being in KCF, and where Uand Vare partitioned conformably with A−λB. We will
now show 1) and 2):
1) Since A−λB is square and of normal rank n−k, it has exactly kleft minimal
indices, say n1, . . . , nkand exactly kright minimal indices, say m1, . . . , mk. We may
assume without loss of generality that they are paired up to form square blocks of
one left and right minimal index each, i.e., we may assume that S(λ) has the block
diagonal form
S(λ) = diag Ln1(λ)>0
0Lm1(λ),...,Lnk(λ)>0
0Lmk(λ).
SOLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION 11
Then the perturbed pencil takes the form
PA−λB +τ(UDAV∗−λUDBV∗)Q=R(λ) 0
τ U2(DA−λDB)V∗
1Rnew(λ),
where
Rnew(λ) := S(λ) + τ U2(DA−λDB)V∗
2.
Clearly, the characteristic polynomial of PA+τ UDAV∗−λ(B+τ UDBV∗)Qis
given by det R(λ)·det Rnew(λ) and from the definition of Rnew(λ) it is clear that the
coefficients of det Rnew(λ) are polynomials in the entries of U2and V∗
2and thus also
of Uand V∗.
Now let λ0be an eigenvalue of R(λ), i.e., det R(λ0) = 0. Note that if ei,j denotes
the ith standard basis vector of Cj, and F= (αi−λ0βi)eni+1,ni+mi+1e∗
ni+1,ni+mi+1,
then
det Lni(λ0)>0
0Lmi(λ0)+F= (−λ0)ni(αi−λ0βi).
Thus, with eithe ith standard basis vector in Cn, for the particular choice
U2=V2= [en1+1 en1+m1+1+n2+1 . . . en1+m1+1+···+nk1+mk−1+1+nk+1]
we obtain that
det Rnew(λ0) = (−λ0)n1+···+nk(α1−λ0β1)· · · (αk−λ0βk)
which is nonzero as the eigenvalues of DA−λDBare by hypothesis distinct from
λ0. But then det Rnew(λ0) is generically nonzero (the set of all (U, V ∗) for which
det Rnew(λ0) = 0 is by definition an algebraic set, because det Rnew(λ0) is a polynomial
in the entries of Uand V∗) which shows that Rnew(λ0) is generically regular and does
not have λ0as an eigenvalue. Since intersections of finitely many generic sets are
still generic we can conclude that the spectra of R(λ) and Rnew(λ) are disjoint. But
then, it immediately follows from [25, Lemma 6.11] and [16, XII.2, Thm. 2] that the
perturbed pencil e
A−λe
Bhas the KCF as given in (4.5).
2) Let λ0be a regular eigenvalue of A−λB and thus of R(λ). It then follows
immediately from 1) that λ0is also an eigenvalue of e
A−λe
Bfor each τ6= 0. For the
moment, let τ6= 0 be fixed and let the columns of Yform a basis of the left null space
Nl(λ0) of e
A−λe
B. Partition
Y∗P−1= [Y∗
1Y∗
2]
conformably with the partition in (4.6). Since λ0is not an eigenvalue of Rnew(λ)
we obtain from Y∗(e
A−λ0e
B) = 0 that Y2= 0. But this implies that the columns
of Yform a basis for the left null space Nl(λ0) for all values of τ∈R\ {0}as the
construction of the transformation matrices Pand Qonly depends on A,B, and U,
but not on τ. Furthermore, we obtain
Y∗U=Y∗P−1PU =Y∗
10·0
U2= 0,
12 HOCHSTENBACH, MEHL, AND PLESTENJAK
i.e., Nl(λ0) is orthogonal to the space spanned by the columns of U.
Observe that the statement on the right null space Nr(λ0) does not follow imme-
diately from the partitioning in (4.6) as in general we have V16= 0. But we can apply
the already proved part of the theorem to the pencil A∗−λB∗and the perturba-
tion V(D∗
A−λD∗
B)U∗to obtain the corresponding statements for the right null space
Nr(λ0). This finishes the proof of 2).
Finally, assume that A−λB does have one of the eigenvalues 0 or ∞. Then apply
a M¨obius transformation of the form
Mα,β(A−λB) := αA +βB −λ(αB −βA)
where α, β ∈Rare such that α2+β2= 1 and such that Mα,β(A−λB) does neither
have the eigenvalues 0 nor ∞. Note that this M¨obius transformation just has the effect
of “rotating” eigenvalues on the extended real line R∪ {∞}, but it leaves eigenvectors
and the Jordan structure invariant, see, e.g., [26]. The result then follows by applying
the already proved parts of the theorem on Mα,β(A−λB) followed by applying the
inverse M¨obius transformation
Mα,−β(C−λD) := αC −βD −λ(αD +βC).
to give the corresponding statements for b
A−λb
B. In particular, this shows 3).
Remark 4.4. We mention that part 1) in Theorem 4.3 is in line with one of the
main results of [7], where it was shown that generically the regular part of a singular
pencil stays invariant under generic perturbations that do not make the pencil regular.
Part 1) of Theorem 4.3 extends this results (in the sense of the theorem) to the case of
rank completing perturbation. Clearly, the regular part of the pencil will be completely
changed if generic perturbations of a rank larger than the difference of the size and
the normal rank of the pencil is applied.
Theorem 4.3 characterizes the properties of the regular eigenvalues of the per-
turbed pencil e
A−λe
Bas in (4.1). We will next investigate the properties of the
eigenvalues from the newly created block Rnew. We start with the following lemma
that will be needed for the main results. The values γ1, . . . , γkin the lemma are the
eigenvalues that we will prescribe later in Theorem 4.6 using the matrices DAand
DB.
Lemma 4.5. Let A−λB be an n×nsingular pencil of normal rank n−kwith
left minimal indices n1, . . . , nkand right minimal indices m1, . . . , mk. Furthermore,
let U, V ∈Cn,k have full column rank, N:= n1+· · ·+nk,M:= m1+· · ·+mk, and let
γ1, . . . , γk∈Cbe given values that are distinct from the regular eigenvalues of A−λB.
Then, generically with respect to the entries of Uand V∗, the following statements
hold:
1. There exist exactly Mpairwise distinct values α1, . . . , αMdifferent from the
eigenvalues of A−λB and different from γ1, . . . , γksuch that for each αithere
exists a nonzero vector ziwith (A−αiB)zi= 0 and V∗zi= 0.
2. There exist exactly Npairwise distinct values β1, . . . , βNdifferent from the
eigenvalues of A−λB and different from γ1, . . . , γkand α1, . . . , αMsuch that
for each βithere exists a nonzero vector wiwith w∗
i(A−βiB)=0and w∗
iUi= 0.
3. For any given set of klinearly independent vectors t1, . . . , tk∈Ckthere exist
nonzero vectors s1, . . . , skwith (A−γiB)si= 0 and ti=V∗sifor i= 1, . . . , k.
SOLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION 13
Proof. 1) Without loss of generality we may assume that A−λB is in KCF such
that the blocks Lm1(λ), . . . , Lmk(λ) associated with the right minimal indices appear
first in the form. Then for each α∈Cdifferent from the eigenvalues of A−λB the
columns of
q1(α). . . qk(α)=
q11(α) 0
...
0qkk(α)
0. . . 0
with qjj(α) =
1
α
.
.
.
αmj
form a basis for ker(A−αB) for each α∈C. (When αis a regular eigenvalue, there
are additional vectors in ker(A−αB) since the rank of A−αB drops below the normal
rank n−k.)
We are looking for z6= 0 and αsuch that V∗z= 0 and (A−αB)z= 0. Since we
want αto be distinct from the regular eigenvalues of A−λB, the vector zhas to be
of the form
z=c1q1(α) + · · · +ckqk(α),
where c= [c1. . . ck]T6= 0. From V∗z= 0 we get the equation
(4.7) G(α)c= 0,
where G(α) is a k×kmatrix whose element gij(α) = v∗
iqj(α) is a polynomial in α
which generically with respect to the entries of v∗
iwill have degree mjfor i, j = 1, . . . , k.
Equation (4.7) has a nontrivial solution if and only if det G(α) = 0, where det G(α) is
a polynomial in αwhich generically with respect to the entries of V∗is of degree M.
Thus det G(α) will have Mroots α1, . . . , αM(counted with multiplicities).
On the other hand, for each fixed µ∈C, we have that det G(µ) is also a polynomial
in the entries of V∗= [v1· · · vk]∗. For the particular choice
v1=e1, v2=em1+2, . . . , vk=em1+···+mk−1+k
we obtain that v∗
iqj(µ) = δij so that G(µ) = Ikshows that det G(µ) is a nonzero poly-
nomial in the entries of V∗. It thus follows that generically with respect to the entries
of V∗we will have det G(µ)6= 0, and consequently the fixed value µwill generically
not be among the roots of G(α) as a polynomial in α. Since the intersection of finitely
many generic sets is still generic, it follows that we can generically exclude finitely
many values from the zeros α1, . . . , αMof G(α). This shows that generically with
respect to the entries of V∗, the values α1, . . . , αMare different from the eigenvalues
of A−λB and also from the given values γ1, . . . , γk.
Next we show that the roots α1, . . . , αMof p(α) := det G(α) generically are pair-
wise distinct which is exactly the case if none of these values is a root of q(α) = p0(α),
where qis the formal derivative of pwith respect to α. Also recall that the polynomials
pand qhave a common root if and only if their resultant is zero, which is given by the
determinant of the Sylvester matrix S(p, q) associated with pand q. Since the entries
of the Sylvester matrix are coefficients of the polynomials pand q, it follows that
det S(p, q) is a polynomial with respect to the entries of V∗. It remains to show that
det S(p, q) is a nonzero polynomial (because then we will have that det S(p, q)6= 0 is
a generic property with respect to the entries of V∗), and for this it is enough to show
14 HOCHSTENBACH, MEHL, AND PLESTENJAK
that for a particular choice of the entries of Vwe have that the values α1, . . . , αM
are pairwise distinct. Now taking v1=em1+1 −ε1e1,v2=em1+m2+2 −ε2em1+2,
. . . , vk=em1+···+mk+k−εkem1+···+mk−1+k, with ε1, . . . , εk>0, we obtain that
v∗
iqj(α) = δijαmj−εjand thus G(α) is diagonal and
det G(α)=(αm1−ε1)· · · (αmk−εk).
Since the roots of each factor (αmj−εj) are mjpairwise distinct complex numbers
lying on a circle centered at zero with radius ε1/mj
j, it remains to choose the values
ε1, . . . , εkin such a way that the kradii are pairwise distinct to obtain the desired
example.
2) In a way similar to the one in 1) we can consider the left null space for A−αB and
show the existence of β1, . . . , βNand the corresponding nonzero vectors w1, . . . , wN,
where now the statements are generic with respect to the entries of U. In particular, by
interpreting Vas already fixed, this shows that generically with respect to the entries
of U, the values β1, . . . , βNare not only different from the eigenvalues of A−λB and
γ1, . . . , γk, but also from the values α1, . . . , αMconstructed in 1).
3) With the same notation as in 1) we now aim to solve the equations
si=c1q1(γi) + · · · +ckqk(γi) and V∗si=ti,
or, equivalently, G(γi)c=tifor i= 1, . . . , k. Since γiis different from the values
α1, . . . , αM, we have det G(γi)6= 0 and hence G(γi)c=tiis uniquely solvable for cfor
i= 1, . . . , k.
The following theorem encapsulates the main result on the new eigenvalues of our
perturbed pencil.
Theorem 4.6. Let A−λB be an n×nsingular pencil of normal rank n−k
with left minimal indices n1, . . . , nkand right minimal indices m1, . . . , mk. Further-
more, let U, V ∈Cn,k have full column rank and let DA=diag(a1, . . . , ak), DB=
diag(b1, . . . , bk)∈Ck,k be such that DA−λDBis regular and such that all (not neces-
sarily pairwise distinct) values γi:= ai
bi,i= 1, . . . , k, are different from the eigenvalues
of A−λB. (Here, ai
biis interpreted as the infinite eigenvalue, if bi= 0.) Finally, let
N:= n1+· · · +nkand M:= m1+· · · +mk. Then generically with respect to the
entries of Uand V∗, the following statements hold:
1. The pencil (4.1) has Msimple eigenvalues α1, . . . , αMwhich are independent
of τ6= 0, so that for each of these eigenvalues its right eigenvector xiis constant
in τ6= 0 (up to scaling) and satisfies V∗xi= 0, while the left eigenvector yiis
a linear function of τ(up to scaling) and satisfies U∗yi6= 0 for all τ6= 0.
2. The pencil (4.1) has Nsimple eigenvalues β1, . . . , βNwhich are independent
of τ6= 0, so that for each of these eigenvalues its left eigenvector yiis constant
for τ6= 0 (up to scaling) and satisfies U∗yi= 0, while the right eigenvector xi
is a linear function of τ(up to scaling) and satisfies V∗xi6= 0 for all τ6= 0.
3. For each τ6= 0 each γiis an eigenvalue of (4.1) with the same algebraic
multiplicity as for the pencil DA−λDB. Furthermore, the left and right null
spaces Nl(γi)and Nr(γi)of (4.1) associated with γiare constant in τ. In
addition, we have:
(a) Nr(γi)∩ker(V∗) = {0}, i.e., for each right eigenvector xof (4.1) associ-
ated with γiwe have V∗x6= 0.
SOLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION 15
(b) Nl(γi)∩ker(U∗) = {0}, i.e., for each left eigenvector yof (4.1) associated
with γiwe have U∗y6= 0.
(Note that the simplicity of the eigenvalues α1, . . . , αM, β1, . . . , βNimplies in particular
that they are all different from the eigenvalues of A−λB and from the eigenvalues
γ1, . . . , γk.)
Proof. Without loss of generality, we may assume that the infinite eigenvalue is not
among the eigenvalues of DA−λDB. Otherwise, we may as in the proof of Theorem 4.3
apply a M¨obius transformation to both A−λB and DA−λDBsuch that DA−λDB
does not have the eigenvalue ∞, apply the statement that was proved for this special
situation, and finally transform back with the inverse M¨obius transformation to obtain
the desired result.
Observe that generically with respect to the entries of Uand V∗, the statements
of Lemma 4.5 will hold, where for the values γ1, . . . , γkwe will take the eigenvalues
of the pencil DA−λDBand for the vectors t1, . . . , tkwe will take the standard basis
vectors e1, . . . , ekfrom Ck. We now show 1)–3).
1) By Lemma 4.5 there exist exactly Mpairwise distinct values α1, . . . , αMdiffer-
ent from the eigenvalues of A−λB and from γ1, . . . , γk, and nonzero vectors z1, . . . , zM
such that (A−αiB)zi= 0 and V∗zi= 0 for i= 1, . . . , M. From this we obtain that
(A−αiB+τ UDAV∗−αiτ UDBV∗)zi= 0 which means that αiis an eigenvalue of
(4.1) for τ6= 0 with a right eigenvector zithat is invariant under τ.
Considering now τas a variable, it follows that the pencil
(4.8) Gi+τHi:= (A−αiB) + τ U(DA−αiDB)V∗
is singular. Suppose that nrank(Gi, Hi) = n−jfor j≥1, which means that (4.8) has j
right and jleft minimal indices. Then we know from Gizi= 0 and Hizi= 0 that one of
these right minimal indices is equal to zero. The remaining j−1 right minimal indices
are all larger than zero, because otherwise there would exist yi∈ker(Gi)∩ker(Hi)
linearly independent of ziwhich implies that αiwould be a multiple eigenvalue of
(4.1) in contradiction to Lemma 4.5.
Now suppose that (4.8) has a left minimal index being zero. Then there exists a
vector wi6= 0 such that w∗
iGi= 0 and w∗
iHi= 0, which implies w∗
iU= 0, because V
was assumed to have full rank and DA−αiDBis nonsingular since the values αiare
different from the eigenvalues of DA−λDB. But then by Lemma 4.5 αiis equal to
one of the values β1, . . . , βNwhich is a contradiction. Thus, all left minimal indices
of (4.8) are larger than or equal to one.
Furthermore, we know that rank(Gi) = n−kand rank(Hi) = ksince αidiffers
from all finite regular eigenvalues of A−λB and all eigenvalues of DA−λDB. It follows
that in the KCF of the pencil (4.8) there are at least n−k−jblocks associated with
the eigenvalue infinity and at least k−jblocks associated with the eigenvalue zero.
By a simple computation we obtain that the dimension of the KCF of (4.8) is at
least (n+j−1) ×(n+j−1); therefore the only option is j= 1 and hence (4.8) has
exactly one right minimal index (being zero) and exactly one left minimal index, say
p. Then another simple computation shows that the dimension of the KCF of (4.8) is
at least (n+p−1) ×(n+p−1) showing that the left minimal index pmust be equal
16 HOCHSTENBACH, MEHL, AND PLESTENJAK
to one. Consequently, there exist linearly independent vectors siand tisuch that
s∗
i(A−αiB) = 0,
t∗
i(A−αiB) + s∗
iU(DA−αiDB)V∗= 0,
t∗
iU(DA−αiDB)V∗= 0
and s∗
iU6= 0. Up to scaling, the left eigenvector yiof (4.1) associated with αithen
has the form yi(τ) = si+τtiand is a linear function of τ.
2) This follows completely analogously to 1).
3) Clearly, the standard basis vectors e1, . . . , ekare eigenvectors of the pencil
DA−λDBassociated with the eigenvalues γ1, . . . , γk. By Lemma 4.5, there exists k
(necessarily linearly independent) vectors s1, . . . , sk∈Cnsuch that (A−γiB)si= 0
and ei=V∗sifor i= 1, . . . , k. Then for each i= 1, . . . , k and each τ6= 0 we have
(e
A−γie
B)si= (A−γiB)si+τ U(DA−γiDB)V∗si= 0.
This implies that the values γ1, . . . , γkare eigenvalues of e
A−λe
Bwith the same al-
gebraic multiplicities as for DA−λDB. Furthermore, it follows that the null space
Nr(γi) does not depend on τand by construction we have Nr(γi)∩ker(V∗) = {0}.
By applying Lemma 4.5 to the pencil A∗−λB∗we obtain the analogous statements
for the left null spaces Nl(γi).
Summary 4.7. Summarizing the results from Theorem 4.3 and Theorem 4.6, let
A−λB be an n×nsingular pencil of normal rank n−k, and with left minimal
indices n1, . . . , nkand right minimal indices m1, . . . , mk, and let U,V,DA,DB,N,
and Mbe as in Theorem 4.6. Since the regular part of A−λB then has the size
r:= n−N−M−kand we have found N+M+knew eigenvalues in Theorem 4.6,
we have classified all eigenvalues of the perturbed pencil
e
A−λe
B:= A−λB +τ(UDAV∗−λ UDBV∗)
into the following three groups:
1. Regular eigenvalues: There are rsuch eigenvalues that are exactly the eigen-
values of A−λB. The corresponding right eigenvectors xand left eigenvectors
ysatisfy V∗x= 0 and U∗y= 0.
2. Prescribed eigenvalues: There are ksuch eigenvalues that coincide with the
keigenvalues of DA−λDB. The corresponding right eigenvectors xand left
eigenvectors ysatisfy both V∗x6= 0 and U∗y6= 0.
3. Random eigenvalues: These are the remaining N+Meigenvalues. They are
simple and if µis such an eigenvalue, xis a corresponding right eigenvector
and yis a corresponding left eigenvector, then we either have V∗x= 0 and
U∗y6= 0, or V∗x6= 0 and U∗y= 0.
Thus, the eigenvalues of A−λB can be identified from the eigenvalues of e
A−λe
B
by investigating orthogonality properties of the corresponding left and right eigenvec-
tors. We will use this observation in the following section for the development of an
algorithm for computing the eigenvalues of a singular square pencil.
Remark 4.8. If Aand Bare symmetric, then it seems that for our current
approach it is necessary to use nonsymmetric rank completing perturbations. Namely,
when a symmetric perturbation of the form τ U(DA−λDB)U∗is used, there is an issue
with the third group in Summary 4.7 as random eigenvalues appear either as double
SOLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION 17
real eigenvalues or in complex conjugate pairs, and in the former case the orthogonality
constraints cannot be satisfied. We leave the study of structured singular pencils for
a future research.
5. A perturbation method for singular generalized eigenvalue problems.
In this section we explain how it is possible in the generic case to extract the finite
regular eigenvalues by solving only one perturbed eigenvalue problem. The key are
the orthogonality properties that are existent or non-existent for the eigenvectors
associated with regular, prescribed, and random eigenvalues, respectively. We can
exploit these properties and obtain reliable results already in double precision.
Let A−λB be a singular n×npencil with a normal rank n−k, where k > 0. We
determine nrank(A, B) by computing rank(A−ζB) for a random ζ. As we have shown
in the previous section, if we take two random n×kmatrices Uand Vwith orthonormal
columns, a regular k×kdiagonal pencil DA−λDB, and τ6= 0, then the perturbed
pencil (4.1) is regular. Regular eigenvalues of A−λB (theoretically) remain constant
under this perturbation. In contrast, eigenvalues that originate from the singular part
of the pencil (the so called “random” eigenvalues) may be “anywhere in the complex
plane”. In addition, (4.1) also has k“prescribed” eigenvalues that coincide with the
eigenvalues of DA−λDB.
Theoretically, if we compute all eigenvalues λitogether with the left and right
eigenvectors xiand yifor i= 1, . . . , n of (4.1), then max(kV∗xik,kU∗yik) = 0 for a
regular eigenvalue and max(kV∗xik,kU∗yik)>0 for a prescribed or a random eigen-
value, so we could apply this criterion to extract the regular eigenvalues. In the
following we will discuss how the above criterion is affected by computation in finite
precision and how does it depend on τ. We will also introduce other criteria that may
be used for the same purpose or to further separate regular eigenvalues into finite and
infinite ones.
If xiand yiare normalized left and right eigenvectors of the perturbed problem
(4.1) for an eigenvalue λi, we can compute the number
(5.1) s(λi) = y∗
ie
Bxi=y∗
iBxi+τ y∗
iUDBV∗xi.
It is easy to see that s(λi)6= 0 for a simple finite eigenvalue λi. As explained in the
following lemma, which is a straightforward generalization of the standard result for a
pencil A−λI, see, e.g., [42, Sec. 2.9], 1/|s(λi)|occurs in the expression for a standard
condition number of a simple finite eigenvalue.
Lemma 5.1. Let λibe a simple finite eigenvalue of a regular matrix pencil e
A−λe
B
and let xiand yibe its normalized left and right eigenvectors. If we perturb the pencil
into (e
A+θE)−λ(e
B+θF)for a small θ > 0, then λiperturbs into
(5.2) λi+θy∗
iExi−λiy∗
iFxi
s(λi)+O(θ2).
For a simple finite regular eigenvalue λiwe have V∗xi= 0 and U∗yi= 0. This
implies s(λi) = y∗
iBxiand consequently s(λi) does not change with τ6= 0. For a
regular infinite eigenvalue we have y∗
iB= 0, Bxi= 0, V∗xi= 0, and U∗yi= 0,
therefore s(∞) = 0, again independent of τ6= 0. On the other hand, we can show
that values s(λ) of prescribed and random eigenvalues depend on τand go to 0 as τ
goes to 0. For this, we need the following lemma.
18 HOCHSTENBACH, MEHL, AND PLESTENJAK
Lemma 5.2. Let A−λB be a singular pencil and let αbe different from all
eigenvalues of A−λB, i.e., rank(A−αB) = nrank(A, B). If y∗(A−αB) = 0 and
(A−αB)x= 0 then y∗Ax =y∗Bx = 0.
Proof. We know from the structure of the left and right singular blocks that
x∈ Nr(α) can be written as x=q1+αq2+· · ·+αpqp+1, where the vectors q1, . . . , qp+1
form the chain
Aq1= 0, Aq2=Bq1, . . . , Aqp+1 =Bqp, Bqp+1 = 0
for certain p≥0. Similarly, y∈ Nl(α) can be written as y=w1+αw2+· · · +αrwr+1,
where the vectors w1, . . . , wr+1 form the chain
A∗w1= 0, A∗w2=B∗w1, . . . , A∗wr+1 =B∗wr, B∗wr+1 = 0
for certain r≥0. To show y∗Bx = 0 it is enough to show that w∗
iBqj= 0 for all
i= 1, . . . , r + 1 and j= 1, . . . , p + 1. For j=p+ 1 this follows from Bqp+1 = 0, so
we can assume that j≤p. It follows that w∗
iBqj=w∗
iAqj+1. If i= 1, then w∗
1A= 0,
if not, we can continue to w∗
iBqj=w∗
iAqj+1 =w∗
i−1Bqj+1.As we continue in this
manner, we eventually reach either w∗
1A= 0 or Bqp+1 = 0. It follows that y∗Bx = 0
and from y∗Ax =αy∗Bx we get y∗Ax = 0 as well.
Lemma 5.3. Let λibe a prescribed or random eigenvalue of (4.1) under the
assumptions of Theorem 4.6. Then there exists a nonzero constant cisuch that s(γi) =
ciτ.
Proof. First, let λibe a prescribed eigenvalue. Then by the proof of Theorem 4.6
the corresponding left and right eigenvectors satisfy y∗
i(A−λiB) = 0 and (A−λiB)xi=
0 which by Lemma 5.2 implies y∗
iBxi= 0. But then we have s(λi) = ciτwith
ci=y∗
iUDBV∗xiand cimust be nonzero, because γiis a simple eigenvalue of e
A−λe
B
for τ6= 0.
Next, let λibe a random eigenvalue, such that V∗xi= 0 and U∗yi6= 0. We know
(see the proof of Theorem 4.6) that yiis a linear function of τas yi(τ) = si+τti,
where s∗
i(A−λiB) = 0. Since (A−λiB)xi= 0, it follows from Lemma 5.2 that
s∗
iBxi= 0 and y∗
ie
Bxi=y∗
iBxi=τ t∗
iBxi. The case V∗xi6= 0 and U∗yi= 0 can be
shown analogously.
So, if we take a small τand if all finite regular eigenvalues are simple and none of
them is too ill-conditioned, then we can separate the finite regular eigenvalues from
the remaining ones with the help of the values s(λ).
Let εbe the machine precision and let the matrices Aand Bbe scaled in such
way that kAk=kBk= 1. If all finite regular eigenvalues are simple and not too
ill-conditioned, then we expect to observe the situation in Table 5.1, where c > 0 is a
constant, independent of τ, and possibly different for each eigenvalue and each entry
in the table.
We now explain the values in the Table 5.1. We will start with column |s(λ)|and a
finite regular eigenvalue, where we assume that all finite regular eigenvalues are simple
and well conditioned. It follows that λis a simple eigenvalue of e
A−λe
B, therefore
y∗e
Bx 6= 0 and, since this value is independent of τand ε, we have |s(λ)|=c. For an
infinite eigenvalue we should have y∗e
Bx = 0 in exact computation, instead, in finite
precision, we get |y∗e
Bx|< cε. Finally, in the generic case, if λiis a prescribed or
SOLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION 19
Table 5.1
Characteristics of the eigenvalues of the perturbed pencil as in (4.1).
Eigenvalue λ|s(λ)| kV∗xk kU∗yk
Finite regular eigenvalue of A−λB c < cε/τ < cε/τ
Infinite regular eigenvalue of A−λB < cε < cε/τ < cε/τ
Prescribed eigenvalue of DA−λDBcτ c c
Random eigenvalue from an Lpblock cτ < cε/τ c
Random eigenvalue from an LT
pblock cτ c < cε/τ
random eigenvalue then it follows from Lemma 5.3 that |s(λi)|=ciτfor a positive
constant ci.
Finally, we turn to the values in the columns kV∗xkand kU∗xkthat are marked
by < cε/τ. In exact arithmetic, these values should be zero. In finite precision
however, due to the supposed backward stability of the applied eigenproblem solver,
the computed eigenvalues and eigenvectors of the matrix pencil e
A−λe
Bare exact
eigenpairs of a perturbed pencil e
A+E−λ(e
B+F), where kEk ≤ c1ke
Akεand kFk ≤
c2ke
Bkε. If we assume that all finite eigenvalues of e
A−λe
Bare simple, then we can say
something about the eigenvector perturbations. This is done in the following lemma
whose proof is omitted, because it is a straightforward generalization of the result for
the pencil A−λI from [42, Sec. 2.10].
Lemma 5.4. Let all finite eigenvalues λiof e
A−λe
Bbe simple and let xiand yi
be corresponding left and right normalized eigenvectors. If the pencil is perturbed into
e
A+θE −λ(e
B+θF), then the eigenvector xiis perturbed into
exi=xi+θ
n
X
k=1,k6=i
y∗
k(E−λiF)xi
(λi−λk)s(λk)xk+O(θ2).
Suppose that λiis a finite regular eigenvalue of A−λB. Then λiis also an
eigenvalue of e
A−λe
Band V∗xi= 0, where xiis an exact normalized right eigenvector.
In finite precision, xibecomes perturbed in directions of other eigenvectors and it
follows from Lemma 5.4 that a contribution in direction of another eigenvector depends
on the condition number of the corresponding eigenvalue. The only contributions
that affect the value of kV∗exikare those related to prescribed eigenvalues or random
eigenvalues from left singular blocks, as right eigenvectors of other eigenvalues are
orthogonal to V. As condition numbers of these eigenvalues are equal to 1/(cτ) and
kV∗xjk=cfor the corresponding right eigenvectors, it follows from Lemma 5.4 and
the backward stability of the computed eigenpairs that kV∗exik< cε/τ.
Next, we now turn to the question what choices for the value τare most beneficial.
If we use a very small τ, i.e., τclose to ε, then prescribed and random eigenvalues are
very ill-conditioned, and perturbations of eigenvectors may affect the values of kV∗xik
and kU∗yikthat much that none of these values may be close to zero, when they
should be. So, if τis too small, we may not be able to use the values of kV∗xikand
kU∗yikto extract the regular eigenvalues. However, if all finite regular eigenvalues are
simple, then we may still use the values |s(λi)|to extract the finite regular eigenvalues.
On the other hand, if we use a large τ, then all eigenvalues (except the infi-
nite one) are expected to be well-conditioned which means that the eigenvectors will
not change much and the computed left and right eigenvectors will be orthogonal to
20 HOCHSTENBACH, MEHL, AND PLESTENJAK
Vor Uin finite precision, when they should be. So, for large τ, we can first use
max(kV∗xik,kU∗yik) to extract the regular eigenvalues and then use |s(λi)|to distin-
guish the finite regular eigenvalues from the infinite one. In practice, we see this as a
better option, because it does not depend on finite regular eigenvalues being simple.
However, we should not choose τtoo large as this may decrease the precision of the
computed finite regular eigenvalues. Since the computed eigenvalues are, due to as-
sumed backward stability, exact eigenvalues of a slightly perturbed pencil e
A−λe
B, it
is safe to use τup to the limit ke
Ak≈kAkand ke
Bk≈kBk.
Based on the above discussion, we summarize our method in Algorithm 1, where
we assume that the pencil is scaled in such way that kAk=kBk= 1. We note that
this scaling is mainly for convenience, to determine an appropriate default value for τ
in Algorithm 1. For efficiency, we can also use the 1-norm scaling kAk1=kBk1= 1.
Algorithm 1: Computing finite regular eigenvalues of a singular pencil
(A, B) by a rank-completing perturbation.
Input: Aand Bsuch that kAk=kBk= 1, perturbation constant τ(default 10−2),
thresholds δ1(default ε1/2) and δ2(default 102ε).
Output: Eigenvalues of the finite regular part.
1: Compute the normal rank of pencil A−λB:k= rank(A−ζB) for random ζ
2: Select random n×kmatrices Uand Vwith orthonormal columns.
3: Select diagonal k×kmatrices DAand DBsuch that the eigenvalues of
(DA, DB) are (likely) different from those of (A, B) (default: choose diagonal
elements of DAand DBuniformly random from the interval [1,2]).
4: Compute the eigenvalues λi,i= 1, . . . , n, and right and left eigenvectors xi
and yiof ( e
A, e
B) = (A+τ UDAV∗, B +τ UDBV∗).
5: Compute si=y∗
ie
Bxifor i= 1, . . . , n.
6: Compute ζi= max(kV∗xik,kU∗yik) for i= 1, . . . , n.
7: Return all eigenvalues λifor i= 1, . . . , n, where ζi< δ1and |si|> δ2.
As we will illustrate by experiments in the next section, the above approach seems
to work very well in double precision for small or moderate singular pencils. Of course,
if some of the eigenvalues are very ill-conditioned (for instance when some of the eigen-
values are multiple), then the method can fail in extracting some of the finite regular
eigenvalues. However, its advantage over staircase-based methods may be the follow-
ing observation: if we make a wrong rank decision in a staircase algorithm, then the
method usually fails completely and returns no eigenvalues at all; see Example 6.4 in
the next section. In contrast, the method proposed here is able to detect (if not all) at
least the well-conditioned finite regular eigenvalues of the pencil under consideration.
6. Numerical examples. In this section we demonstrate the method with sev-
eral numerical examples computed in Matlab 2015b.
Example 6.1. For the first example we take a 7 ×7 pencil A−λB, where
A=
−1−1−1−1−1−1−1
1000000
1211111
1233333
1232222
1234333
1234554
, B =
−2−2−2−2−2−2−2
2−1−1−1−1−1−1
2555555
2554444
2556555
2556777
2556766
.
SOLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION 21
The matrices are built in such way that the KCF of the pencil contains blocks of all
four possible types, nrank(A, B) = 6 and the pencil is singular. Its KCF has blocks
J1(1/2), J2(1/3), N1,L1, and LT
2. If we apply Algorithm 1, we get the values in the
following table:
k λk|sk| kV∗xkk kU∗ykk
1 0.333333 1.5·10−21.3·10−15 1.3·10−14
2 0.500000 9.5·10−41.3·10−14 1.9·10−14
3∞3.8·10−19 2.8·10−15 1.3·10−14
4−0.244794 + 0.421723i7.8·10−35.8·10−25.6·10−15
5−0.244794 −0.421723i7.8·10−35.8·10−25.6·10−15
6 0.383682 2.1·10−42.6·10−24.2·10−1
7 0.478292 2.6·10−49.2·10−15 5.2·10−1
The values in the table follow the pattern from the previous section and it is easy to
detect that λ1and λ2are finite regular eigenvalues, λ3is a regular infinite eigenvalue,
λ4, λ5, and λ7are random eigenvalues, and λ6is the prescribed eigenvalue.
Example 6.2. For the next example we take example C3 from [12]. This example
comes from control theory and belongs to a set of examples C1, C2, and C3, where
each has successively more ill-conditioned eigenvalues. The pencil has the form
A−λB =
1−2 100 0 0
1 0 −1 0 0
0 0 0 1 −75
0 0 0 0 2
−λ
01000
00100
00010
00001
.
Its KCF contains blocks L2,J1(1), and J1(2). As the pencil is rectangular, we add a
zero line to make it square. This adds an LT
0block to the KCF. Algorithm 1 returns
the following table for A−λB, from which the finite regular eigenvalues 1 and 2 can
be extracted.
k λk|sk| kV∗xkk kU∗ykk
1 1.000000 1.2·10−21.7·10−15 1.9·10−15
2 2.000000 1.2·10−22.3·10−15 1.7·10−15
3−0.693767 + 1.563033i2.1·10−22.3·10−15 5.0·10−1
4−0.693767 −1.563033i2.1·10−22.3·10−15 5.0·10−1
5 78.673901 2.8·10−33.3·10−16.4·10−1
As in [12] we add some noise and perturb initial A−λB into b
A−λb
Bby adding
10−6rand(4,5) to Aand B. Regular eigenvalues of b
A−λb
Bcan still be extracted by
Algorithm 1 if we adjust the parameter δ1. The values we get are in the following
table:
k λk|sk| kV∗xkk kU∗ykk
1 0.999990 7.6·10−32.6·10−15 5.2·10−7
2 2.000058 7.6·10−32.8·10−15 5.6·10−7
3 101.850555 8.2·10−41.3·10−13.7·10−1
4−14.308508 1.9·10−25.8·10−16 4.5·10−1
5 15.734162 9.3·10−31.0·10−17 4.7·10−1
Example 6.3. This is an example from [14, Sec. 5], where the staircase algorithm
fails to find a regular subspace of proper size under a small random perturbation. We
take
A1−λB1=
0010
0001
0000
−λ
δ000
0δ0 0
0100
,
where δ= 1.5·10−8. The KCF structure of the pencil is J2(0) and L1which means
that 0 is a double regular eigenvalue. It is reported in [14] that if we add a random
22 HOCHSTENBACH, MEHL, AND PLESTENJAK
perturbation of size 10−14 to the pencil, then Guptri reports the regular part J1(0) and
we could confirm this using a Matlab implementation of Guptri in [34]. If we enlarge
the perturbation to 10−11, Guptri returns no regular part at all, while Algorithm 1
returns two finite regular eigenvalues λ1and λ2from the following table.
k λk|sk| kV∗xkk kU∗ykk
1−1.4306543 ·10−31.6·10−11 5.5·10−17 6.7·10−10
29.9599790 ·10−41.6·10−11 0.0 6.7·10−10
3−2.2641370 ·1075.2·10−92.9·10−18 2.6·10−6
4 1.1878888 ·1001.6·10−32.0·10−17.8·10−1
Example 6.4. For this example we take the singular pencil ∆1−λ∆0of size
300 ×300 from [33, Ex. 18]. This example is related to two random matrices Aand
Bof size 10 ×10 in a way that the finite regular eigenvalues of ∆1−λ∆0are exactly
the values λsuch that A+λB has a multiple eigenvalue (see Section 3.2). We know
from the properties of the problem, that there are 90 such values λand that the KCF
of ∆1−λ∆0contains 100 N1and 10 left and 10 right singular blocks. The conjecture
from [33] is that the singular blocks are 5 LT
4, 5 LT
5, 5 L5, and 5 L6blocks.
This example is also available as demo_double_eig_mp in toolbox MultiParEig
[35]. The staircase algorithm in MultiParEig fails to extract the finite regular part
of size 90 in double precision, but manages to extract all 90 finite regular eigenvalues
using quadruple precision and the Multiprecision Computing Toolbox [36]. If we apply
Algorithm 1 on ∆1−λ∆0in double precision, we get the following values:
k λk|sk| kV∗xkk kU∗ykk
1 0.508999 + 2.016378i3.0·10−32.3·10−14 1.6·10−14
.
.
..
.
..
.
..
.
..
.
.
89 4.266290 −0.925962i1.4·10−62.6·10−13 7.2·10−14
90 −0.628208 3.2·10−72.8·10−14 1.3·10−11
91 ∞1.1·10−17 7.1·10−15 7.1·10−15
.
.
..
.
..
.
..
.
..
.
.
190 ∞2.8·10−21 5.9·10−15 7.9·10−15
191 −6.276934 3.2·10−72.7·10−14 4.5·10−5
.
.
..
.
..
.
..
.
..
.
.
300 7.125982 2.3·10−51.7·10−11.0·10−2
From the columns kV∗xkkand kU∗ykkwe get maxk=1,...,190(max(kV∗xkk,kU∗ykk)) =
1.3·10−11 and mink=191,...,300 max(kV∗xkk,kU∗ykk)=4.5·10−5, which shows a clear
gap that separates regular eigenvalues from the prescribed and random ones. Next,
in the set of regular eigenvalues there is also a clear gap between s90 and s91 that
separates finite regular eigenvalues from infinite ones, since mink=1,...,90 |sk|= 3.2·10−7
and maxk=91,...,190 |sk|= 1.1·10−17.
7. The singular two-parameter eigenvalue problem. In a two-parameter
eigenvalue problem (2EP) [1] we have the equations
(A1+λB1+µC1)x1= 0,
(7.1) (A2+λB2+µC2)x2= 0,
where A1,B1, and C1are of size n1×n1, and A2,B2, and C2are of size n2×n2.
Sought are scalars λ, µ and nonzero vectors x1and x2such that (7.1) is satisfied.
We say that (λ, µ) is an eigenvalue of the 2EP and the tensor product x1⊗x2is the
SOLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION 23
corresponding eigenvector. Define the operator determinants
∆0=B1⊗C2−C1⊗B2,
∆1=C1⊗A2−A1⊗C2,(7.2)
∆2=A1⊗B2−B1⊗A2.
Then problem (7.1) is related to a coupled pair of GEPs
∆1z=λ∆0z,
(7.3) ∆2z=µ∆0z
for a decomposable tensor z=x1⊗x2. If ∆0is nonsingular, then Atkinson [1] shows
that the solutions of (7.1) and (7.3) agree and the matrices ∆−1
0∆1and ∆−1
0∆2com-
mute. In the nonsingular case the two-parameter problem (7.1) has n1n2eigenvalues
and it can be solved with a variant of the QZ algorithm on (7.3); see [18].
It turns out that for many problems occurring in practice both pencils (∆1,∆0)
and (∆2,∆0) are singular and we have a singular two-parameter eigenvalue problem
[31]. Applications involve delay-differential equations [20], quadratic two-parameter
eigenvalue problems [32, 19], model updating [5], roots of systems of bivariate poly-
nomials [37, 2], and others.
There exists a staircase type algorithm that works on both singular pencils (7.3)
simultaneously and extracts regular finite eigenvalues; see [32] and an implementation
in [35]. However, as we have shown in Examples 6.3 and 6.4, a staircase algorithm
can fail. In this section we propose another method that can be applied to a singular
2EP, which in some cases finds regular finite eigenvalues when the staircase algorithm
fails, while in some other cases the situation is exactly the opposite.
We can apply Algorithm 1 to ∆1z=λ∆0z, one out of two singular pencils in
(7.3) to compute λicomponents of eigenvalues (λi, µi). This is, however, only half
of the required information and for each λiwe have to find the corresponding µi.
Subsequently, we insert λ=λiin (7.1) and search for common eigenvalues µof a pair
of pencils (A1−λiB1)−µC1and (A2−λiB2)−µC2that could be singular as well. We
detect the common eigenvalues by comparing the sets of computed eigenvalues for the
first and the second pencil, for which we use Algorithm 1 again. The overall method
is given in Algorithm 2.
Algorithm 2: Computing finite regular eigenvalues of a singular 2EP
Input: Matrices A1, B1, C1, A2, B2, C2from (7.1) which also provide ∆1and ∆0from
(7.2); threshold δ(default δ=ε1/2), and parameters for Algorithm 1.
Output: Finite regular eigenvalues and eigenvectors of (7.1).
1: Compute finite regular eigenvalues λ1, . . . , λrof ∆1−λ∆0using Alg. 1.
2: for j= 1, . . . , r
3: Compute eigenvalues µ(1)
1, . . . , µ(1)
m1of (A1−λjB1)−µC1using Alg. 1.
4: Compute eigenvalues µ(2)
1, . . . , µ(2)
m2of (A2−λjB2)−µC2using Alg. 1.
5: Reorder eigenpairs: |µ(1)
1−µ(2)
1| ≤ · · · ≤ |µ(1)
m−µ(2)
m|for m= min(m1, m2).
6: for k= 1, . . . , m
7: if |µ(1)
k−µ(2)
k|< δ then add (λj,1
2(µ(1)
k+µ(2)
k)) to list of eigenvalues.
Some remarks:
24 HOCHSTENBACH, MEHL, AND PLESTENJAK
•If we know that each eigenvalue has a unique λcomponent, then we can replace
Lines 6 and 7 by selecting (λj,1
2(µ(1)
1+µ(2)
1)), x(1)
1⊗x(2)
1) regardless of the
difference |µ(1)
1−µ(2)
1|.
•If n1=n2=nthen the complexity of Line 1 is O(n6) while the complexity of
Lines 2 to 7 is at most O(n5) in case r=O(n2).
Example 7.1. Consider a system of bivariate polynomials
p1(λ, µ) = 1 + 2λ+ 3λ+ 4λ2+ 5λµ + 6µ2+ 7λ3+ 8λ2µ+ 9λµ2+ 10µ3= 0,
p2(λ, µ) = 10 + 9λ+ 8µ+ 7λ2+ 6λµ + 5µ2+ 4λ3+ 3λ2µ+ 2λµ2+µ3= 0.
Using a uniform determinantal representation from [2], we write the above system as
a two parameter eigenvalue problem of the form
A1+λB1+µC1=
0 0 4 + 7λ1 0
0 5 + 8λ2−λ1
6+9λ+ 10µ3 1 0 −λ
1−µ0 0 0
0 1 −µ0 0
,
A2+λB2+µC2=
0 0 7 + 4λ1 0
0 6 + 3λ9−λ1
5+2λ+µ8 10 0 −λ
1−µ0 0 0
0 1 −µ0 0
,
where pi(λ, µ) = det(Ai+λBi+µCi) for i= 1,2. The obtained two-parameter
eigenvalue problem is singular and has 9 regular eigenvalues (λj, µj) which are exactly
the 9 solutions of the initial polynomial system.
If we apply Algorithm 2 to the above problem, we get all 9 solutions. In Line 2 we
compute first components λ1, . . . , λ9as finite regular eigenvalues of the corresponding
singular pencil ∆1−λ∆0from (7.3), whose KCF contains 4 L0, 4 LT
0, 2 N4, 1 N2, 2
N1, and 9 J1blocks. For each λjwe compute the candidates for µjin Lines 4 and 5,
where the KCF of singular pencils (Ai−λjBi)−µCicontains 1 N2and 3 J1blocks
for i= 1,2 and j= 1,...,9.
We remark that the above approach might also fail, in particular if we apply
it to systems of bivariate polynomials of high degree. It can happen that some of
the eigenvalues of ∆1−λ∆0are so ill-conditioned that the algorithm cannot separate
them from the infinite eigenvalues. In such a case a possible solution would be to apply
computation in higher precision, using, e.g., the Multiprecision Computing Toolbox
[36].
8. Conclusions. We have proposed a method to approximate the finite eigen-
values of a singular pencil by means of a so-called rank-completing perturbation, i.e.,
a random perturbation of the pencil of rank k, where n−kis the normal rank of the
pencil of dimension n×n. The use of such perturbations ensures that, generically, the
regular finite and infinite eigenvalues remain fixed, while there appear newly generated
eigenvalues. For most problems we can well distinguish the original regular eigenval-
ues from the newly created ones by looking at the angles of the eigenvectors with
respect to the perturbation spaces, and at the condition numbers of the eigenvalues.
Therefore, this method may be useful for a wide range of applications.
When the eigenvalues of the original problem are very ill-conditioned, it may
become difficult to distinguish between original eigenvalues and newly created ones,
and it is an outcome to have another alternative type of methods available: the class
SOLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION 25
of staircase algorithms, such as e.g., Guptri [17] or a staircase type algorithm for
singular two-parameter eigenvalue problems [32] in [35]. These methods can be rapid
and accurate for many problems. However, the key part of staircase techniques are
a number of rank decisions, which can be difficult and ill-posed, see e.g., [14] and
Examples 6.3 and 6.4. In such cases, these methods may not return even a single
eigenvalue.
A code for the approach developed in this paper can be obtained via [35].
Acknowledgment: The authors would like to thank Stefan Johansson for providing a
beta version of Matrix Canonical Structure (MCS) Toolbox [34] that includes a Matlab
implementation of Guptri and is an important alternative to the original Guptri [17]
that we can no longer use in Matlab due to the 32-bit limitation.
Genealogical acknowledgment: During this research project the first two authors
found out that they are twelfth cousins. Christian and Michiel thank their common
ancestors Caspar H¨olterhoff (1552–1625) and Catharina Teschemacher (ca. 1555–1639)
for making this possible.
REFERENCES
[1] F. V. Atkinson,Multiparameter Eigenvalue Problems, Academic Press, New York, 1972.
[2] A. Boralevi, J. van Doornmalen, J. Draisma, M. E. Hochstenbach, and B. Plesten-
jak,Uniform determinantal representations, SIAM J. Appl. Algebra Geometry, 1 (2017),
pp. 415–441.
[3] P. Benner, P. Losse, V. Mehrmann, and M. Voigt,Numerical Linear Algebra Methods
for Linear Differential-Algebraic Equations. In: A. Ilchmann, T. Reis (eds), Surveys in
Differential-Elgebraic Equations III. Differ.-Algebr. Equa. Forum, pp. 117-175, Springer,
Cham, 2015.
[4] R. Byers, C. He, and V. Mehrmann,Where is the nearest non-regular pencil?, Linear Algebra
Appl., 285 (1998), pp. 81–105.
[5] N. Cottin,Dynamic model updating—a multiparameter eigenvalue problem, Mechanical systems
and signal processing, 15 (2001), pp. 649–665.
[6] E.J. Davison and S.H. Wang,Properties and calculation of transmission zeros of linear mul-
tivariable systems, Automatica, 10 (1974), pp. 643-658.
[7] F. De Ter´
an and F. M. Dopico,Low rank perturbation of Kronecker structures without full
rank, SIAM J. Matrix Anal. Appl., 29 (2007), pp. 496–529.
[8] , A note on generic Kronecker orbits of matrix pencils with fixed rank, SIAM J. Matrix
Anal. Appl., 30 (2008), pp. 491–496.
[9] F. De Ter´
an, F. M. Dopico, and J. Moro,First order spectral perturbation theory of square
singular matrix pencils, Linear Algebra Appl., 429 (2008), pp. 548–576.
[10] J. Demmel,Generalized Non-Hermitian Eigenproblems. Section 2.6 in: Z. Bai, J. Demmel,
J. Dongarra, A. Ruhe, and H. van der Vorst (eds), Templates for the Solution of Algebraic
Eigenvalue Problems: a practical guide, pp. 28–36, SIAM, Philadelphia, 2000.
[11] J. Demmel and B. K˚
agstr¨
om,Computing stable eigendecompositions of matrix pencils, Linear
Algebra Appl., 88/89 (1987), pp. 139–186.
[12] , Accurate solutions of ill-posed problems in control theory, SIAM J. Matrix. Anal. Appl.,
9 (1988), pp. 126–145.
[13] , The generalized Schur decomposition of an arbitrary pencil A−λB–robust software with
error bounds and applications. Part I: theory and algorithms, ACM Trans. Math. Software,
19 (1993), pp. 160–174.
[14] A. Edelman and Y. Ma,Staircase failures explained by orthogonal versal forms, SIAM J. Ma-
trix Anal. Appl., 21 (2000), pp. 1004–1025.
[15] A. Emami-Naeini and P. Van Dooren,Computation of zeros of linear multivariable systems,
Automatica 18 (1982), pp. 415-430.
[16] F.R. Gantmacher,Theory of Matrices. Volume 1 and 2, Chelsea, New York, 1959.
[17] Guptri, software for singular pencils, www8.cs.umu.se/research/nla/singular_pairs/guptri/.
[18] M. E. Hochstenbach, T. Koˇ
sir, and B. Plestenjak,A Jacobi–Davidson type method for
26 HOCHSTENBACH, MEHL, AND PLESTENJAK
the nonsingular two-parameter eigenvalue problem, SIAM J. Matrix Anal. Appl., 26 (2005),
pp. 477–497.
[19] M. E. Hochstenbach, A. Muhiˇ
c, and B. Plestenjak,On linearizations of the quadratic
two-parameter eigenvalue problem, Linear Algebra Appl., 436 (2012), pp. 2725–2743.
[20] E. Jarlebring and M. E. Hochstenbach,Polynomial two-parameter eigenvalue problems
and matrix pencil methods for stability of delay-differential equations, Linear Algebra Appl.,
431 (2009), pp. 369–380.
[21] E. Jarlebring, S. Kvaal, and W. Michiels,Computing all pairs (λ,µ) such that λis a double
eigenvalue of A+µB, SIAM J. Matrix Anal. Appl., 32 (2011), pp. 902–927.
[22] B. K˚
agstr¨
om,Singular matrix pencils. Section 8.7 in: Z. Bai, J. Demmel, J. Dongarra, A. Ruhe,
and H. van der Vorst (eds), Templates for the Solution of Algebraic Eigenvalue Problems: a
practical guide, pp. 260–277, SIAM, Philadelphia, 2000.
[23] P. Kunkel and V. Mehrmann,Differential-Algebraic Equations: Analysis and Numerical
Solution, EMS Publishing House, Z¨urich, Switzerland, 2006.
[24] A.J. Laub and B.C. Moore,Calculation of transmission zeros using QZ techniques, Auto-
matica, 14 (1978), pp. 557–566.
[25] D. S. Mackey, N. Mackey, C. Mehl, and V. Mehrmann,Skew-symmetric matrix polyno-
mials and their Smith forms, Linear Algebra Appl., 438 (2013), pp. 4625–4653.
[26] , M¨obius transformations of matrix polynomials, Linear Algebra Appl., 470 (2015), pp. 120–
184.
[27] C. Mehl, V. Mehrmann, and M. Wojtylak,On the distance to singularity via low rank
perturbations, Operators and Matrices, 9 (2015), pp. 733–772.
[28] , Parameter-dependent rank-one perturbations of singular Hermitian or symmetric pencils,
SIAM J. Matrix Anal. Appl., 38 (2017), pp. 72–95.
[29] , Linear algebra properties of dissipative Hamiltonian descriptor systems, Preprint 1-2018,
Institut f¨ur Mathematik, TU Berlin (2018), submitted.
[30] V. Mehrmann, private communication, 2018.
[31] A. Muhiˇ
c and B. Plestenjak,On the singular two-parameter eigenvalue problem, Elec-
tron. J. Linear Algebra, 18 (2009), pp. 420–437.
[32] , On the quadratic two-parameter eigenvalue problem and its linearization, Linear Algebra
Appl., 432 (2010), pp. 2529–2542.
[33] , A method for computing all values λsuch that A+λB has a multiple eigenvalue, Linear
Algebra Appl., 440 (2014), pp. 345–359.
[34] MCS Toolbox. The Matrix Canonical Structure Toolbox for Matlab, www.cs.umu.se/english/
research/groups/matrix-computations/stratigraph.
[35] MultiParEig. Toolbox for multiparameter eigenvalue problems, www.mathworks.com/
matlabcentral/fileexchange/47844-multipareig.
[36] Multiprecision Computing Toolbox. Advanpix, Tokyo. www.advanpix.com.
[37] B. Plestenjak and M. E. Hochstenbach,Roots of bivariate polynomial systems via deter-
minantal representations, SIAM J. Sci. Comput., 38 (2016), pp. A765–A788.
[38] N. Valeev,On a spectral property of irregular pencils, Ufa Mathematical Journal, 4 (2012),
pp. 44–52.
[39] , On quasiregular spectrum of matrix pencils, Doklady Mathematics, vol. 88, Springer,
2013, pp. 545–547.
[40] P. Van Dooren,The computation of Kronecker’s canonical form of a singular pencil, Linear
Algebra Appl., 27 (1979), pp. 103–140.
[41] , Reducing subspaces: Definitions, properties and algorithms. In: B. K˚agstr¨om, A. Ruhe
(eds), Matrix Pencils. Lecture Notes in Mathematics, vol 973. Springer, Berlin, Heidelberg,
1983, pp. 58–73.
[42] J.H. Wilkinson,Algebraic Eigenvalue Problem, Oxford University Press, Oxford, 1965.
[43] , Kronecker’s canonical form and the QZ algorithm, Linear Algebra Appl., 28 (1979),
pp. 285–303.
