Solving singular generalized eigenvalue problems by a rank-completing perturbation
aa r X i v : . [ m a t h . NA ] A p r SOLVING SINGULAR GENERALIZED EIGENVALUE PROBLEMSBY 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 withrespect to accuracy and efficiency. The existing package Guptri is very elegant but may be time-demanding, evenfor 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 isboth 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 perturbation, Guptri,model updating, double eigenvalues, two-parameter eigenvalue problem, differential algebraic equations, quadratictwo-parameter eigenvalue problem.
AMS subject classifications.
1. Introduction.
We study the computation of eigenvalues of small to medium-sized matrixpencils A − λB , where A and B are (real or complex) n × m matrices such that the matrix pencil A − λB is singular, which means that m = n , or if m = n thendet( A − λB ) ≡ . In these cases, the common definition of eigenvalues as roots of det( A − λB ) would only bemeaningful for the case m = n , but turns out to be useless as any value λ ∈ C would be aneigenvalue. Therefore, finite eigenvalues of a singular matrix pencil A − λB are typically definedas values λ ∈ C satisfying rank( A − λ B ) < nrank( A, B ), wherenrank(
A, B ) := max ζ ∈ C rank( A − ζB )denotes the normal rank of the pencil A − λB ; see [10]. Similarly, we say that ∞ is an eigenvalueof the singular pencil A − λB if rank( B ) < nrank( A, B ). In the following we will mainly restrictourselves to the case m = n as the case m = n can easily be reduced to the square case by addingan appropriate number of zero rows or columns.The singular generalized eigenvalue problem (singular GEP) is well known to be ill-conditionedas arbitrarily small perturbation may cause drastic changes in the eigenvalues. A classical exampleis given by the pencils A − λB = (cid:20) (cid:21) − λ (cid:20) (cid:21) and e A − λ e B = (cid:20) ε ε (cid:21) − λ (cid:20) ε ε (cid:21) , where ε , . . . , ε ∈ C \ { } ; see [22]. While A − λB is singular and has only the eigenvalue 1, theperturbed pencil e A − λ e B is regular and has the eigenvalues ε ε and ε ε that can be anywhere inthe complex plane even for tiny absolute values of ε , . . . , ε . ∗ Version April 3, 2019. Department of Mathematics and Computer Science, TU Eindhoven, PO Box 513, 5600MB, The Netherlands, ∼ hochsten . This author has been supported by an NWO Vidi researchgrant. † Institut f¨ur Mathematik, Technische Universit¨at Berlin, Sekretariat MA 4-5, Straße des 17. Juni 136, 10623Berlin, Germany, [email protected] . This author has been supported by a Dutch 4TU AMI visitor’s grant. ‡ IMFM and Faculty of Mathematics and Physics, University of Ljubljana, Jadranska 19, SI-1000 Ljubljana,Slovenia, [email protected] . This author has been supported in part by the Slovenian ResearchAgency (grant P1-0294). 1
HOCHSTENBACH, MEHL, AND PLESTENJAK
On the other hand, it was observed in [43] that a situation as above is exceptional and thatgenerically small perturbations of a singular square pencil make the pencil regular and some ofthe eigenvalues of the perturbed pencil are very close to the original eigenvalues of the singularpencil. 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 U and V are certain random 6 × We see that the three (finite) eigenvalues of the regular part are correct. Following the terminologyof [40], the other three values are “fake eigenvalues” and correspond to the singular part of thepencil. (Explicit error analysis for the eigenvalues of singular pencils has been undertaken in[9, 11].) Despite this observation, Van Dooren suggests in [40] to solve the singular generalizedeigenvalue problem by first extracting the regular part and then use the QZ algorithm on thatpart. Wilkinson strongly supports this recommendation in [43].A robust software package which follows Van Dooren’s approach is Guptri [13, 17]. For asingular pencil, 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 usually excellent, this method may be quite time-consuming; forinstance, applying Guptri on a singular 300 ×
300 pencil on our machine took over 20 seconds,while Matlab’s eig on a random pencil of the same size spent less than a second. Another issueis the fact that staircase type methods such as Guptri need rank decisions. If the pencil has aminimal index of size η (see Section 2 for more details), then at least η + 1 such decisions have tobe taken. Typically, these decisions tend to become more and more critical during a run of thestaircase algorithm. See, e.g., [14] or [33, Ex. 18], where a variant of the staircase algorithm forthe singular two-parameter eigenvalue problem, introduced in [31], fails in double precision butgets the right result in higher precision.Another way of extracting the regular part using fewer rank decisions has been suggested in[30]. One may view the singular pencil as a constant coefficient differential-algebraic equationand perform a regularization procedure with the help of a derivative array as described in [3]. Inthis way, the regular part of the pencil can be extracted by 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 thus result in highcomputational costs.We propose a new method to compute the eigenvalues of a singular pencil. The method isbased on considering perturbations of rank k = n − nrank( A, B )which we will call rank-completing perturbations as the rank is exactly large enough to genericallyturn the pencil into a pencil of full normal rank. As we will show, the canonical form of the originalregular part of the given pencil stays invariant under generic rank-completing perturbations.The idea of computing eigenvalues of singular pencils with rank-completing perturbationsis not completely new, and the following specific type has been used in system theory as earlyas in the 70s (without the use of the terminology “rank-completing perturbation”). If a linear We note that this experiment has been performed some years ago. A current practical issue is that there is nopublicly available 64-bit Guptri code.OLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION x = Ax + Bu,y = Cx + Du is given, where A ∈ R n,n , B ∈ R n,m , C ∈ R r,n , and D ∈ R r,m are the system matrices, x standsfor the state of the system, u is the input, and y is the output, then the eigenvalues of the systempencil S ( λ ) = (cid:20) λI − A B − C D (cid:21) are of particular interest in control theory; see [15] and the references therein. (If the system isminimal, then these eigenvalues are also referred to as transmission zeros of the system.) Clearly,if m = r , then the pencil S ( λ ) is rectangular and thus singular. For that case and under theadditional assumptions r < m and nrank S ( λ ) = n + r , the following algorithm based on ideas of[6] has been proposed in [24] for the computation of the transmission zeros: Select random matrices [ E F ], [ E F ] ∈ R m − r,n + m so that S i ( λ ) := λI − A B − C DE i F i is regular for i = 1 , Compute the eigenvalues E i of S i ( λ ) for i = 1 , Compute the intersection E = E ∩ E .Since for each eigenvalue λ of S ( λ ) we have rank S i ( λ ) < n + m , it immediately followsthat the eigenvalues of S ( λ ) are contained in the spectrum of S i ( λ ) for both i = 1 ,
2. Theextended matrix pencils will give rise to two sets of fake eigenvalues. As generically these setswill be disjoint if the applied perturbations are generated randomly, it follows that the set E willgenerically coincide with the set of eigenvalues of S ( λ ).However, as pointed out in [15], this method may encounter difficulties in distinguishingthe finite zeros from the infinite ones, in particular if the latter occur with a high multiplicity.Another problem may occur in identifying the values that belong to the intersection E . Althoughthe original eigenvalues of the pencil theoretically coincide with a subset of both E and E , theymay still differ slightly in practice due to finite precision arithmetic. Therefore, a tolerance has tobe prescribed that decides when two values are considered to be equal. If this tolerance is chosentoo small, then some of the eigenvalues may be missed. If, on the other hand, the tolerance isset too large, then two close fake eigenvalues of S ( λ ) and S ( λ ) may be falsely identified as aneigenvalue of S ( λ ).In this paper, we show that the eigenvalues of a singular pencil can be efficiently computedwith the help of just one rank-completing perturbation of the form A − λB + τ U ( D A − λD B ) V ∗ , where U , V are n × k matrices with orthonormal columns, D A , D B are diagonal k × k matrices, and τ is a nonzero scalar. The orthonormality of the columns is not strictly necessary, but convenient,for instance since in this case the norm of the perturbation can easily be controlled by theparameter τ . The problem of identifying the subset of eigenvalues of the original pencil among thecomputed eigenvalues of the perturbed pencil is then taken care of by the key observation that theleft and right eigenvectors that correspond to the true eigenvalues satisfy orthogonality relations HOCHSTENBACH, MEHL, AND PLESTENJAK with respect to the matrices U and V . Thus, instead of comparing the spectra of two differentpencils, the true eigenvalues can be separated from the fake eigenvalues by using informationfrom the corresponding left and right eigenvectors from only one perturbed pencil. We note thatperturbations of singular matrix pencils have already been considered in [4, 9, 27, 38, 39], butit 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 −
2. Preliminaries.
We will interpret matrix pencils both as pairs of matrices (
A, B ) ∈ C n,m × C n,m or as n × m matrix polynomials A − λB of degree at most one and we will switch betweenthese notations whenever useful. An important tool in the theory of singular pencils is theKronecker 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 × m matrix pencil.Then there exist nonsingular matrices P ∈ C n,n and Q ∈ C m,m such that (2.1) P ( A − λB ) Q = (cid:20) R ( λ ) 00 S ( λ ) (cid:21) , R ( λ ) = (cid:20) J − λI r I s − λN (cid:21) with J , N in Jordan normal form and in addition N being nilpotent, and S ( λ ) = diag (cid:0) L m ( λ ) , . . . , L m k ( λ ) , L n ( λ ) ⊤ , . . . , L n ℓ ( λ ) ⊤ (cid:1) , where L j ( λ ) = [0 I j ] − λ [ I j is of size j × ( j + 1) , and m i ≥ for i = 1 , . . . , k , and n i ≥ for 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 J are exactly the finite eigenvalues of A − λB , while theeigenvalue 0 of N corresponds to the infinite eigenvalue of A − λB . The parameters m , . . . , m k and n , . . . , n ℓ are called the right and the left minimal indices of A − λB , respectively. One mayeasily check that the normal rank nrank( A, B ) is equal to min( n − ℓ, m − k ). In the remainderof the paper, we will consider square pencils A − λB , i.e., we have n = m . Note that this implies k = ℓ , i.e., we must have the same number of right and left minimal indices. However, theparticular values of the left and right minimal indices may be distinct.In contrast to the eigenvalues of singular pencils, the corresponding eigenvectors and deflatingsubspaces are not well defined. To understand why, we consider the following example borrowedfrom [29] (and slightly adapted). The pencil(2.2) A − λB = − λ − λ
10 0 0 obviously has the regular part R ( λ ) = [1 − λ ] and thus the pencil A − λB has the single eigenvalue λ = 1 with algebraic multiplicity one. Nevertheless, any vector of the form x ( α, β ) := [ α β β ] ⊤ with α, β ∈ C satisfies Ax = λ Bx and thus could be interpreted as an eigenvector of the pencil.One may argue that the choice α = 0 and β = 0 seems to be canonical and gives a unique OLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION /α − β/α − λ − λ
10 0 0 α β β = − λ − λ
10 0 0 that for any choice of α, β with α = 0, the vector x ( α, β ) can be used to extract the regular partof the pencil as well (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 donot consider corresponding eigenvectors. However, as we will see in Section 4, the eigenvectors ofthe perturbed pencil play a key role in our approach.Instead of eigenvectors and deflating subspaces, the concept of reducing subspaces introducedin [41] is more adequate in the case of singular pencils. We say that a subspace M is a reducingsubspace for the pencil A − λB if dim( A M + B M ) = dim( M ) − k , where k is the number of rightsingular blocks. In the example above, the reducing subspace associated with the eigenvalue λ = 1 is exactly given by all vectors x ( α, β ) with α, β ∈ C . The minimal reducing subspace M RS ( A, B ) is the intersection of all reducing subspaces. It is unique and can be numericallycomputed 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 = A r − λB r × × A reg − λB reg × A l − λB l , where the matrices P and Q are unitary, A r − λB r has only right singular blocks in its KCF, A l − λB l has only left singular blocks in its KCF, and A reg − λB reg has only regular blocks in itsKCF. We will briefly come back to a minimal reducing subspace in Section 3.3.
3. Motivation and applications.
Before proposing our new method, we first review somemotivating applications where one is interested in computing the eigenvalues of a singular pencil,besides the computation of transmission zeros already mentioned in the introduction.
Linear differential al-gebraic equations (DAEs) with constant coefficients have the general form E ˙ x = Ax + f ( t ) , x ( t ) = x , where E, A ∈ R k,n , t ∈ R , x ∈ R n , and f : [ t , ∞ ) → R k is a given inhomogeneity; see [23].Linear time-invariant descriptor systems consist of a DAE combined with a system input andoutput and take the form E ˙ x = Ax + Bu, x ( t ) = x ,y = Cx + Du, where, in addition, B ∈ R k,m , C ∈ R p,n , D ∈ R p,m . Here, x : [ t , ∞ ) → R n stands for the stateof the descriptor system, u : [ t , ∞ ) → R m is the input and y : [ t , ∞ ) → R p the output. Ashighlighted in [3], the problem may be well-posed even if the underlying pencil A − λE is singular.Indeed, even in the singular case the corresponding DAE may have a solution (even unique) forparticular inhomogeneities f and special initial conditions ( t , x ). HOCHSTENBACH, MEHL, AND PLESTENJAK
Given two n × n matrices A and B , we are interested inall values λ such that A + λB has a double eigenvalue. For the generic case of a double eigenvalue,we look for independent vectors x and y such that( A + λB − µI ) x = 0 , ( A + λB − µI ) y = 0 , for some µ . 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 a linearization of the second equation and followed by thesolution of the obtained singular two-parameter eigenvalue problem. The wanted values λ areeigenvalues of the singular pencil ∆ − λ ∆ of size 3 n × n , where∆ = A ⊗ R − I ⊗ P, ∆ = B ⊗ R − I ⊗ Q, for P = A AB + BA − A I
00 0 I , Q = B − B − I , R = − B I − I . Generically, the pencil ∆ − λ ∆ will have normal rank 3 n − n . For more details, see [33] as wellas [19] for other possible linearizations. In the singular two-parameter eigen-value problem we are given a pair of singular pencils ∆ − λ ∆ and ∆ − µ ∆ (see also (7.2) and(7.3)), and the goal is to find finite regular eigenvalues ( λ , µ ), where (see, e.g., [31] for moredetails): • λ is a finite eigenvalue of ∆ − λ ∆ , • µ is a finite eigenvalue of ∆ − µ ∆ , • and there exists a nonzero vector z such that (∆ − λ ∆ ) z = 0, (∆ − µ ∆ ) z = 0, and z
6∈ M RS (∆ i , ∆ ) for i = 1 ,
4. Rank-completing perturbations of singular pencils.
Let A − λB be a singularpencil, where A, B ∈ C n,n and nrank( A, B ) = n − k . In this section we investigate the effect of rank-completing perturbations , i.e., rank- k generic perturbations of the form(4.1) e A − λ e B := A − λB + τ ( U D A V ∗ − λ U D B V ∗ ) , where D A , D B ∈ C k,k are diagonal matrices such that D A − λD B is a regular pencil, U, V ∈ C n,k have full column rank, and τ ∈ C is nonzero. We investigate the dependence of eigenvalues andeigenvectors of the perturbed pencils on τ .Above and in the following, the term generic is understood in the following sense. A set A ⊆ C m is called algebraic if it is the set of common zeros of finitely many complex polynomialsin m variables, and A is called proper if A 6 = C m . A set Ω ⊆ C m is called generic if its complementis contained in a proper algebraic set. In this sense we say that a property P holds genericallywith respect to the entries of U and V , if there exists a generic set Ω ⊆ ( C n,k ) (where we interpret( C n,k ) as C nk ) such that P holds 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 perturbations of rank k . Indeed, completing U and V to nonsingular matrices P = [ U e U ] and Q = [ V e V ], we obtain OLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION P − ( U D A V ∗ − λ U D B V ∗ )( Q ∗ ) − = (cid:20) D A − λD B
00 0 (cid:21) , which means that the perturbation pencil U D A V ∗ − λ U D B V ∗ has the regular part D A − λD B of size k × k , while, generically, a matrix pencil of nrank k < n would have no regular part [8].Thus, more generally one could consider perturbations of the form(4.2) ( A + U V ∗ + V U ∗ , B + U W ∗ + W U ∗ ) , where U , V , W ∈ C n,ℓ and U , V , W ∈ C n,k − ℓ have full column rank for ℓ ∈ { , , . . . , k } .However, we will restrict ourselves to perturbations of the form (4.1), because of their favorableproperties.Generically, a rank completing perturbation (4.1) of A − λB will result in a regular perturbedpencil e A − λ e B . We will show in the following, that generically the KCF of e A − λ e B is for all τ = 0given by R reg ( λ ) 0 00 R pre ( λ ) 00 0 R ran ( λ ) , where R reg ( λ ) is the regular part of the original pencil A − λB , R pre ( λ ) = D A − λD B and R ran ( λ )only has simple eigenvalues that are different from the eigenvalues of R reg ( λ ) and R pre ( λ ). Thus,the eigenvalues of e A − λ e B are exactly the (say) p eigenvalues of the original pencil A − λB (counted with multiplicities) and n − p newly generated eigenvalues which consist of k “prescribed”eigenvalues which are the eigenvalues of the perturbation pencil U ( D A − λD B ) V ∗ , and n − p − k “random” eigenvalues that are gathered in the part R ran ( λ ).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 willstate some results in more generality covering also the case of rectangular pencils. The followingproposition is a generalization of [28, Thm. 4.2] (which deals with Hermitian pencils) to a blockcase without any specific structure in A and B . Proposition 4.2.
Let A − λB be an n × m singular matrix pencil having at least k leftminimal indices, and let U ∈ C n,k . Then generically (with respect to the entries of U ) there existnonsingular matrices P, Q such that P ( A − λB ) Q = (cid:20) R ( λ ) 00 S ( λ ) (cid:21) and P U = (cid:20) e U (cid:21) , where R ( λ ) and S ( λ ) are the regular and singular parts of A − λB , respectively, and P U ispartitioned conformably with P ( A − λB ) Q .Proof . Without loss of generality we may assume that A − λB is already in the KCF A − λB = (cid:20) R ( λ ) 00 S ( λ ) (cid:21) , R ( λ ) = (cid:20) J − λI I − λN (cid:21) , where S ( λ ) is the singular part and R ( λ ) is the regular part of A − λB with J, N in Jordannormal canonical form and, in addition, N nilpotent. Our main strategy is to use the part of P that corresponds to k arbitrarily chosen left minimal indices to introduce zeros in the componentsof U that correspond to the regular part of A − λB . Here, we can treat the components of U HOCHSTENBACH, MEHL, AND PLESTENJAK corresponding 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 eigenvaluepart” is completely analogous. Since we will only transform parts of the singular part S ( λ )that correspond to k arbitrarily 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 k left 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 U and A − λB have the forms U = kr U n +1 U ... ... n k +1 U k and A − λB = r n . . . n k r J − λI n +1 L n ( λ ) ⊤ ... . . . n k +1 L n k ( λ ) ⊤ . In the following, we will use the notation G ℓ = [0 I ℓ ] ⊤ and H ℓ = [ I ℓ ⊤ to write L ℓ ( λ ) ⊤ = G ℓ − λH ℓ for a singular block corresponding to a left minimal index ℓ . For the transformationmatrices P, Q we make the ansatz P = r n +1 . . . n k +1 r I P . . . P kn +1 I ... . . . n k +1 I and Q = r n . . . n k r I − P H n . . . − P k H n k n I ... . . . n k I As P and Q are designed in such a way that B remains unchanged, we get P ( A − λB ) Q = J − λI P G n − J P H n . . . P k G n k − J P k H n k L n ( λ ) ⊤ . . . L n k ( λ ) ⊤ . It remains to find solutions P i to the equations P i G n i − J P i H n i = 0 for i = 1 , . . . , k to obtain P ( A − λB ) Q = A − λB . Setting P i = [ p i, p i, . . . p i,n i ] for i = 1 , . . . , k , the equations to besolved take the form[ p i, . . . p i,n i ] = P i G n i = J P i H n i = [ J p i, . . . J p i,n i − ] , i = 1 , . . . , k. Thus, we may choose P i = [ p i, J p i, . . . J n i p i, ], where p i, ∈ C r is arbitrary. We will nowuse the freedom in the choice of p i, to guarantee that P U has the desired form. To this end, let u ⊤ i, , . . . , u ⊤ i,n i ∈ C k be the rows of U i , i = 1 , . . . , k . Then the first block component of P U is givenby U + P U + · · · + P k U k and to make it zero, we have to solve the equation(4.3) − U = k X i =1 P i U i = k X i =1 n i X j =0 J j p i, u ⊤ i,j for p , , . . . , p k, . Using the vec-operation that “vectorizes” a matrix by stacking the columnson top of each other and recalling the well-known identity vec( XY Z ) = ( Z ⊤ ⊗ X ) vec( Y ) formatrices X, Y, Z , where ⊗ denotes the Kronecker product, we obtain that(4.4) − vec( U ) = k X i =1 n i X j =0 (cid:0) u i,j ⊗ J j (cid:1) p i, = M · (cid:20) p , ... p k, (cid:21) , OLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION M = h n X j =0 u ,j ⊗ J j . . . n k X j =0 u k,j ⊗ J j i ∈ C rk,rk . The determinant of M is a polynomial in the nk entries of U (in fact, it only depends on theentries of U , . . . , U k ) which is nonzero for the particular choice u , = e , . . . , u k, = e k and u i,j = 0 for j >
0, where e , . . . , e k denote the standard basis vectors of C k . (Indeed, in thiscase M is just the identity of size rk × rk .) Thus, generically (with respect to the entries of U ), the matrix M is invertible, so equation (4.4) and thus also (4.3) can be uniquely solved for p , , . . . , p k, which finishes the proof.The next result shows that the canonical form of the regular part of the original pencil staysinvariant under a generic rank-completing perturbation of the form (4.1). Concerning the eigen-values of the perturbed pencil that are also eigenvalues of the original singular pencil, the resultalso states that the corresponding left and right eigenvectors satisfy a particular orthogonalityrelation. Theorem 4.3.
Let A − λB be an n × n singular pencil of normal rank n − k , let U, V ∈ C n,k have full column rank and let D A , D B ∈ C k,k be such that D A − λD B is regular and all eigenvaluesof D A − λD B are distinct from the eigenvalues of A − λB . Then, generically with respect to theentries of U and V ∗ , the following statements hold for the pencil (4.1) :1. For each τ = 0 , there exist nonsingular matrices e P and e Q such that (4.5) e P (cid:0) e A − λ e B (cid:1) e Q = (cid:20) R ( λ ) 00 R new ( λ ) (cid:21) , where R ( λ ) is the regular part of the original pencil A − λB , and R new ( λ ) is regular andall its eigenvalues are distinct from the eigenvalues of R ( λ ) .2. If λ is a finite eigenvalue of A − λB , i.e., rank ( A − λ B ) < n − k , then λ is an eigenvalueof (4.1) for each τ = 0 . Furthermore, the right null space N r ( λ ) := ker (cid:0) e A − λ e B (cid:1) andthe left null space N l ( λ ) := ker (cid:0)(cid:0) e A − λ e B (cid:1) ∗ (cid:1) are both constant in τ = 0 . In addition:(a) N r ( λ ) ⊥ span( V ) , i.e., if x is a right eigenvector of e A − λ e B associated with λ , then V ∗ x = 0 .(b) N l ( λ ) ⊥ span( U ) , i.e., if y is a left eigenvector of e A − λ e B associated with λ , then U ∗ y = 0 .3. If ∞ is a eigenvalue of A − λB , i.e., rank ( B ) < n − k , then ∞ is an eigenvalue of (4.1) for each τ = 0 . The right and left null spaces N r ( ∞ ) := ker( e B ) and N l ( ∞ ) := ker (cid:0) e B ∗ (cid:1) are both constant in τ = 0 . In addition:(a) N r ( ∞ ) ⊥ span( V ) , i.e., if x is a right eigenvector of e A − λ e B associated with ∞ , then V ∗ x = 0 .(b) N l ( ∞ ) ⊥ span( U ) , i.e., if y is a left eigenvector of e A − λ e B associated with ∞ , then U ∗ y = 0 .Proof . First, we will assume that A − λB does not have one of the eigenvalues 0 or ∞ andwe will show 1) and 2) for this particular case.Applying Proposition 4.2, there generically exist nonsingular matrices P, Q such that(4.6) P ( A − λB ) Q = (cid:20) R ( λ ) 00 S ( λ ) (cid:21) , P U = (cid:20) U (cid:21) , Q ∗ V = (cid:20) V V (cid:21) , where R ( λ ) and S ( λ ) are the regular and singular parts of A − λB , respectively, both being inKCF, and where U and V are partitioned conformably with A − λB . We will now show 1) and2):0 HOCHSTENBACH, MEHL, AND PLESTENJAK
1) Since A − λB is square and of normal rank n − k , it has exactly k left minimal indices,say n , . . . , n k and exactly k right minimal indices, say m , . . . , m k . We may assume without lossof generality that they are paired up to form square blocks of one left and right minimal indexeach, i.e., we may assume that S ( λ ) has the block diagonal form S ( λ ) = diag (cid:18)(cid:20) L n ( λ ) ⊤ L m ( λ ) (cid:21) , . . . , (cid:20) L n k ( λ ) ⊤ L m k ( λ ) (cid:21)(cid:19) . Then the perturbed pencil takes the form P ( e A − λ e B ) Q = (cid:20) R ( λ ) 0 τ U ( D A − λD B ) V ∗ R new ( λ ) (cid:21) , where R new ( λ ) := S ( λ ) + τ U ( D A − λD B ) V ∗ . Clearly, the determinant of P ( e A − λ e B ) Q is equal to det R ( λ ) · det R new ( λ ) and from the definitionof R new ( λ ) it is clear that the coefficients of det R new ( λ ) are polynomials in the entries of U and V ∗ and thus also of U and V ∗ .Now let λ be an eigenvalue of R ( λ ), i.e., det R ( λ ) = 0. Note that if e j,ℓ denotes the j thstandard basis vector of C ℓ , and F i = ( α i − λ β i ) e n i +1 ,n i + m i +1 e ∗ n i +1 ,n i + m i +1 , thendet (cid:18)(cid:20) L n i ( λ ) ⊤ L m i ( λ ) (cid:21) + F i (cid:19) = ( − λ ) n i ( α i − λ β i ) . Thus, with e j the j th standard basis vector in C n , for the particular choice U = V = [ e n +1 e n + m +1+ n +1 . . . e n + m +1+ ··· + n k − + m k − +1+ n k +1 ]we obtain thatdet R new ( λ ) = τ k ( − λ ) n + ··· + n k ( α − λ β ) · · · ( α k − λ β k )which is nonzero as the eigenvalues of D A − λD B are by hypothesis distinct from λ . Butthen det R new ( λ ) is generically nonzero (the set of all ( U, V ∗ ) for which det R new ( λ ) = 0 is bydefinition an algebraic set, because det R new ( λ ) is a polynomial in the entries of U and V ∗ )which shows that R new ( λ ) is generically regular and does not have λ as an eigenvalue. Sinceintersections of finitely many generic sets are still generic we can conclude that the spectra of R ( λ ) and R new ( λ ) 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 B has the KCF as given in (4.5).2) Let λ be a eigenvalue of A − λB and thus of R ( λ ). It then follows directly from 1) that λ is also an eigenvalue of e A − λ e B for each τ = 0. For the moment, let τ = 0 be fixed and let thecolumns of Y form a basis of the left null space N l ( λ ) of e A − λ e B . Partition Y ∗ P − = [ Y ∗ Y ∗ ]conformably with the partition in (4.6). Since λ is not an eigenvalue of R new ( λ ) we obtain from Y ∗ ( e A − λ e B ) = 0 that Y = 0. But this implies that the columns of Y form a basis for the leftnull space N l ( λ ) for all values τ = 0 as the construction of the transformation matrices P and Q only depends on A , B , and U , but not on τ . Furthermore, we obtain Y ∗ U = Y ∗ P − P U = [ Y ∗ · (cid:20) U (cid:21) = 0 , OLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION N l ( λ ) is orthogonal to the space spanned by the columns of U .Observe that the statement on the right null space N r ( λ ) does not follow immediately fromthe partitioning in (4.6) as in general we have V = 0. But we can apply the already provedpart of the theorem to the pencil A ∗ − λB ∗ and the perturbation V ( D ∗ A − λD ∗ B ) U ∗ to obtain thecorresponding statements for the right null space N r ( λ ). This finishes the proof of 2).Finally, assume that A − λB does have one of the eigenvalues 0 or ∞ . Then apply a M¨obiustransformation of the formM α,β ( A − λB ) := αA + βB − λ ( αB − βA )where α, β ∈ R are such that α + β = 1 and such that M α,β ( A − λB ) does neither have theeigenvalues 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 structureinvariant, see, e.g., [26]. The result then follows by applying the already proved parts of thetheorem on M α,β ( A − λB ) followed by applying the inverse M¨obius transformationM α, − β ( C − λD ) := αC − βD − λ ( αD + βC ) . to give the corresponding statements for A − λ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 resultsof [7], where it was shown that generically the regular part of a singular pencil stays invariantunder generic perturbations that do not make the pencil regular. Part 1) of Theorem 4.3 extendsthis result (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 largerthan the difference of the size and the normal rank of the pencil are applied.Theorem 4.3 characterizes the properties of the eigenvalues from the block R of the perturbedpencil e A − λ e B as in (4.1), i.e., of the eigenvalues that coincide with the eigenvalues of theunperturbed pencil. We will next investigate the properties of the eigenvalues from the newlycreated block R new . We start with the following lemma that will be needed for the main results.The values γ , . . . , γ k in the lemma are the eigenvalues that we will prescribe later in Theorem 4.6using the matrices D A and D B . Lemma 4.5.
Let A − λB be an n × n singular pencil of normal rank n − k with left minimalindices n , . . . , n k and right minimal indices m , . . . , m k . Furthermore, let U, V ∈ C n,k have fullcolumn rank, N := n + · · · + n k , M := m + · · · + m k , and let γ , . . . , γ k ∈ C be given valuesthat are distinct from the eigenvalues of A − λB . Then, generically with respect to the entries of U and V ∗ , the following statements hold:1. There exist exactly M pairwise distinct values α , . . . , α M different from the eigenvaluesof A − λB and different from γ , . . . , γ k such that for each α i there exists a nonzero vector z i with ( A − α i B ) z i = 0 and V ∗ z i = 0 .2. There exist exactly N pairwise distinct values β , . . . , β N different from the eigenvalues of A − λB and different from γ , . . . , γ k and α , . . . , α M such that for each β i there exists anonzero vector w i with w ∗ i ( A − β i B ) = 0 and w ∗ i U = 0 .3. For any given set of k linearly independent vectors t , . . . , t k ∈ C k there exist nonzerovectors s , . . . , s k with ( A − γ i B ) s i = 0 and t i = V ∗ s i for i = 1 , . . . , k .Proof . 1) Without loss of generality we may assume that A − λB is in KCF such that theblocks L m ( λ ) , . . . , L m k ( λ ) associated with the right minimal indices appear first in the form.2 HOCHSTENBACH, MEHL, AND PLESTENJAK
Then for each α ∈ C different from the eigenvalues of A − λB the columns of (cid:2) q ( α ) . . . q k ( α ) (cid:3) = q ( α ) 0. . .0 q kk ( α )0 . . . with q jj ( α ) = α ... α m j form a basis for ker( A − αB ). (When α is an eigenvalue of A − λB , there are additional vectorsin ker( A − αB ) since the rank of A − αB drops below the normal rank n − k .)We are looking for z = 0 and α such that V ∗ z = 0 and ( A − αB ) z = 0. Since we want α tobe distinct from the eigenvalues of A − λB , the vector z has to be of the form z = c q ( α ) + · · · + c k q k ( α ) , where c = [ c . . . c k ] T = 0. From V ∗ z = 0 we get the equation(4.7) G ( α ) c = 0 , where G ( α ) is a k × k matrix whose element g ij ( α ) = v ∗ i q j ( α ) is a polynomial in α which genericallywith respect to the entries of v ∗ i will have degree m j for i, j = 1 , . . . , k . Equation (4.7) hasa nontrivial solution if and only if det G ( α ) = 0, where det G ( α ) is a polynomial in α whichgenerically with respect to the entries of V ∗ is of degree M . Thus det G ( α ) will have M roots α , . . . , α M (counted with multiplicities).On the other hand, for each fixed µ ∈ C , we have that det G ( µ ) is also a polynomial in theentries of V ∗ = [ v . . . v k ] ∗ . For the particular choice v = e , v = e m +2 , . . . , v k = e m + ··· + m k − + k we obtain that v ∗ i q j ( µ ) = δ ij so that G ( µ ) = I k shows that det G ( µ ) is a nonzero polynomial inthe entries of V ∗ . It thus follows that generically with respect to the entries of V ∗ we will havedet G ( µ ) = 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 followsthat we can generically exclude finitely many values from the zeros α , . . . , α M of G ( α ). Thisshows that generically with respect to the entries of V ∗ , the values α , . . . , α M are different fromthe eigenvalues of A − λB and also from the given values γ , . . . , γ k .Next we show that the roots α , . . . , α M of p ( α ) := det G ( α ) generically are pairwise distinct.This is exactly the case if the discriminant Disc( p ) of p is nonzero. Since Disc( p ) is a polynomialin the entries of p (this is well known, but can also be seen from the fact that the discriminantis a scalar multiple of the determinant of the Sylvester matrix S ( p, p ′ ) associated with p and itsformal derivative p ′ ), it follows that Disc( p ) is a polynomial with respect to the entries of V ∗ . Itremains to show that Disc( p ) is a nonzero polynomial (because then we will have that Disc( p ) = 0is a generic property with respect to the entries of V ∗ ), and for this it is enough to show that for aparticular choice of the entries of V we have that the values α , . . . , α M are pairwise distinct. Nowtaking v = e m +1 − ε e , v = e m + m +2 − ε e m +2 , . . . , v k = e m + ··· + m k + k − ε k e m + ··· + m k − + k ,with ε , . . . , ε k >
0, we obtain that v ∗ i q j ( α ) = δ ij α m j − ε j and thus G ( α ) is diagonal anddet G ( α ) = ( α m − ε ) · · · ( α m k − ε k ) . Since the roots of each factor ( α m j − ε j ) are m j pairwise distinct complex numbers lying on acircle centered at zero with radius ε /m j j , it remains to choose the values ε , . . . , ε k in such a waythat the k radii are pairwise distinct to obtain the desired example. OLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION A − αB andshow the existence of β , . . . , β N and the corresponding nonzero vectors w , . . . , w N , where nowthe statements are generic with respect to the entries of U . In particular, by interpreting V asalready fixed, this shows that generically with respect to the entries of U , the values β , . . . , β N are not only different from the eigenvalues of A − λB and γ , . . . , γ k , but also from the values α , . . . , α M constructed in 1).3) With the same notation as in 1) we now aim to solve the equations s i = c q ( γ i ) + · · · + c k q k ( γ i ) and V ∗ s i = t i , or, equivalently, G ( γ i ) c = t i for i = 1 , . . . , k . Since γ i is different from the values α , . . . , α M , wehave det G ( γ i ) = 0 and hence G ( γ i ) c = t i is uniquely solvable for c for i = 1 , . . . , k .The following theorem encapsulates the main result on the new eigenvalues of our perturbedpencil. Theorem 4.6.
Let A − λB be an n × n singular pencil of normal rank n − k with left minimalindices n , . . . , n k and right minimal indices m , . . . , m k . Furthermore, let U, V ∈ C n,k have fullcolumn rank and let D A = diag ( a , . . . , a k ) , D B = diag ( b , . . . , b k ) ∈ C k,k be such that D A − λD B is regular and such that all (not necessarily pairwise distinct) values γ i := a i b i , i = 1 , . . . , k , aredifferent from the eigenvalues of A − λB . (Here, a i b i is interpreted as the infinite eigenvalue, if b i = 0 .) Finally, let N := n + · · · + n k and M := m + · · · + m k . Then generically with respectto the entries of U and V ∗ , the following statements hold:1. The pencil (4.1) has M simple eigenvalues α , . . . , α M which are independent of τ = 0 ,so that for each of these eigenvalues its right eigenvector x i is constant in τ = 0 (up toscaling) and satisfies V ∗ x i = 0 , while the left eigenvector y i is a linear function of τ (upto scaling) and satisfies U ∗ y i = 0 for all τ = 0 .2. The pencil (4.1) has N simple eigenvalues β , . . . , β N which are independent of τ = 0 ,so that for each of these eigenvalues its left eigenvector y i is constant for τ = 0 (up toscaling) and satisfies U ∗ y i = 0 , while the right eigenvector x i is a linear function of τ (upto scaling) and satisfies V ∗ x i = 0 for all τ = 0 .3. For each τ = 0 each γ i is an eigenvalue of (4.1) with the same algebraic multiplicity asfor the pencil D A − λD B . Furthermore, the left and right null spaces N l ( γ i ) and N r ( γ i ) of (4.1) associated with γ i are constant in τ . In addition, we have:(a) N r ( γ i ) ∩ ker( V ∗ ) = { } , i.e., for each right eigenvector x of (4.1) associated with γ i we have V ∗ x = 0 .(b) N l ( γ i ) ∩ ker( U ∗ ) = { } , i.e., for each left eigenvector y of (4.1) associated with γ i we have U ∗ y = 0 .(Note that the simplicity of the eigenvalues α , . . . , α M , β , . . . , β N implies that they are all dif-ferent from the eigenvalues of A − λB and γ , . . . , γ k .)Proof . Without loss of generality, we may assume that the infinite eigenvalue is not amongthe eigenvalues of D A − λD B . Otherwise, we may as in the proof of Theorem 4.3 apply a M¨obiustransformation to both A − λB and D A − λD B such that D A − λD B does not have the eigenvalue ∞ , apply the statement that was proved for this special situation, and finally transform backwith the inverse M¨obius transformation to obtain the desired result.Observe that generically with respect to the entries of U and V ∗ , the statements of Lemma 4.5hold, if we take the eigenvalues of the pencil D A − λD B for the values γ , . . . , γ k and take standardbasis vectors e , . . . , e k from C k for the vectors t , . . . , t k . We now show 1)–3).1) By Lemma 4.5 there exist exactly M pairwise distinct values α , . . . , α M different from theeigenvalues of A − λB and from γ , . . . , γ k , and nonzero vectors z , . . . , z M such that ( A − α i B ) z i =0 and V ∗ z i = 0 for i = 1 , . . . , M . From this we obtain that ( A − α i B + τ U D A V ∗ − α i τ U D B V ∗ ) z i =4 HOCHSTENBACH, MEHL, AND PLESTENJAK α i is an eigenvalue of (4.1) for τ = 0 with a right eigenvector z i that is invariantunder τ .Considering now τ as a variable, it follows that the pencil(4.8) G i + τ H i := ( A − α i B ) + τ U ( D A − α i D B ) V ∗ is singular. Suppose that nrank( G i , H i ) = n − j for j ≥
1, which means that (4.8) has j rightand j left minimal indices. We know from G i z i = 0 and H i z i = 0 that one right minimal indexis equal to zero. The remaining j − y i ∈ ker( G i ) ∩ ker( H i ) linearly independent of z i which implies that α i would 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 w i = 0 such that w ∗ i G i = 0 and w ∗ i H i = 0, which implies w ∗ i U = 0, because V has full rank and D A − α i D B is nonsingular since α i is different from the eigenvalues of D A − λD B . But then byLemma 4.5 α i is equal to one of the values β , . . . , β N which is a contradiction. Thus, all leftminimal indices of (4.8) are larger than or equal to one.Furthermore, we know that rank( G i ) = n − k and rank( H i ) = k since α i differs from all finiteeigenvalues of A − λB and all eigenvalues of D A − λD B . It follows that in the KCF of the pencil(4.8) there are at least n − k − j blocks associated with the eigenvalue infinity and at least k − j blocks 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 − × ( n + j − j = 1 and hence (4.8) has exactly one right minimalindex (being zero) and exactly one left minimal index, say p . Then another simple computationshows that the dimension of the KCF of (4.8) is at least ( n + p − × ( n + p −
1) showing thatthe left minimal index p must be equal to one. Consequently, there exist linearly independentvectors w i and z i such that w ∗ i ( A − α i B ) = 0 ,z ∗ i ( A − α i B ) + w ∗ i U ( D A − α i D B ) V ∗ = 0 ,z ∗ i U ( D A − α i D B ) V ∗ = 0and w ∗ i U = 0. Up to scaling, the left eigenvector y i of (4.1) associated with α i then has the form y i ( τ ) = w i + τ z i and is a linear function of τ .2) This follows completely analogously to 1).3) Clearly, the standard basis vectors e , . . . , e k are eigenvectors of the pencil D A − λD B associated with the eigenvalues γ , . . . , γ k . By Lemma 4.5, there exist k (necessarily linearlyindependent) vectors s , . . . , s k ∈ C n such that ( A − γ i B ) s i = 0 and e i = V ∗ s i for i = 1 , . . . , k .Then we have( e A − γ i e B ) s i = ( A − γ i B ) s i + τ U ( D A − γ i D B ) V ∗ s i = 0for each τ = 0 for i = 1 , . . . , k . This implies that the values γ , . . . , γ k are eigenvalues of e A − λ e B with the same algebraic multiplicities as for D A − λD B . Furthermore, it follows that the nullspace N r ( γ i ) does not depend on τ and by construction we have N r ( γ i ) ∩ ker( V ∗ ) = { } .By applying Lemma 4.5 to the pencil A ∗ − λB ∗ we obtain the analogous statements for theleft null spaces N l ( γ i ). Summary 4.7.
Summarizing the results from Theorem 4.3 and Theorem 4.6, let A − λB be an n × n singular pencil of normal rank n − k with left minimal indices n , . . . , n k and rightminimal indices m , . . . , m k , and let U , V , D A , D B , N , and M be as in Theorem 4.6. Since the OLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION A − λB then has size r := n − N − M − k and we have found N + M + k neweigenvalues in Theorem 4.6, we have classified all eigenvalues of the perturbed pencil e A − λ e B := A − λB + τ ( U D A V ∗ − λ U D B V ∗ )into the following three groups:1. True eigenvalues : There are r such eigenvalues that are exactly the eigenvalues of A − λB .The corresponding right eigenvectors x and left eigenvectors y satisfy V ∗ x = 0 and U ∗ y =0.2. Prescribed eigenvalues : There are k such eigenvalues that coincide with the k eigenvaluesof D A − λD B . The corresponding right eigenvectors x and left eigenvectors y satisfy both V ∗ x = 0 and U ∗ y = 0.3. Random eigenvalues : These are the remaining N + M eigenvalues. They are simple andif µ is such an eigenvalue with the corresponding right eigenvector x and left eigenvector y , then we either have V ∗ x = 0 and U ∗ y = 0, or V ∗ x = 0 and U ∗ y = 0.Thus, the eigenvalues of A − λB can be identified from the eigenvalues of e A − λ e B by investigatingorthogonality properties of the corresponding left and right eigenvectors. We will use this obser-vation in the following section for the development of an algorithm for computing the eigenvaluesof a singular square pencil. Remark 4.8. If A and B are symmetric, then it seems that for our current approach we haveto use nonsymmetric rank completing perturbations. Namely, when a symmetric perturbationof the form τ U ( D A − λD B ) U ∗ is used, there is an issue with the third group in Summary 4.7as random eigenvalues appear either as double real eigenvalues or in complex conjugate pairs,and in the former case the orthogonality constraints cannot be satisfied. We leave the study ofstructured singular pencils for future research.
5. A perturbation method for singular generalized eigenvalue problems.
In thissection we explain how in the generic case we can extract the finite true eigenvalues numericallyeven in double precision by solving only one perturbed eigenvalue problem. The key is formed bythe existent or non-existent orthogonality properties of the left and right eigenvectors associatedwith true, prescribed, and random eigenvalues, respectively.Let A − λB be a singular n × n pencil with normal rank n − k , where k >
0. We determinenrank(
A, B ) by computing rank( A − ζB ) for a random ζ . As we have shown in the previoussection, if we take two random n × k matrices U and V with orthonormal columns, a regular k × k diagonal pencil D A − λD B , and τ = 0, then the perturbed pencil (4.1) is regular. The“true” eigenvalues of A − λB (theoretically) remain constant under this perturbation. In contrast,eigenvalues that originate from the singular part of the pencil (the “random” eigenvalues) maybe “anywhere in the complex plane”. In addition, (4.1) also has k “prescribed” eigenvalues thatcoincide with the eigenvalues of D A − λD B .In theory, if we compute all eigenvalues λ i together with the left and right eigenvectors x i and y i for i = 1 , . . . , n of (4.1), then max( k V ∗ x i k , k U ∗ y i k ) = 0 for a true eigenvalue andmax( k V ∗ x i k , k U ∗ y i k ) > τ . We will also introduce othercriteria that may be used for the same purpose or to further separate true eigenvalues into finiteand infinite ones.If x i and y i are normalized left and right eigenvectors of the perturbed problem (4.1) for aneigenvalue λ i , we can compute the number(5.1) s ( λ i ) = y ∗ i e Bx i = y ∗ i Bx i + τ y ∗ i U D B V ∗ x i . HOCHSTENBACH, MEHL, AND PLESTENJAK
It is easy to see that s ( λ i ) = 0 for a simple finite eigenvalue λ i . As explained in the followinglemma, 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 simplefinite eigenvalue. Lemma 5.1.
Let λ i be a simple finite eigenvalue of a regular matrix pencil e A − λ e B and let x i and y i be its normalized left and right eigenvectors. If we perturb the pencil into ( e A + θE ) − λ ( e B + θF ) for a small θ > , then λ i perturbs into (5.2) λ i + θ y ∗ i Ex i − λ i y ∗ i F x i s ( λ i ) + O ( θ ) . If λ i is a simple finite true eigenvalue, then V ∗ x i = 0 and U ∗ y i = 0, which implies that s ( λ i ) = y ∗ i Bx i does not change with τ = 0. For a regular infinite eigenvalue we have y ∗ i B = 0, Bx i = 0, V ∗ x i = 0, and U ∗ y i = 0, therefore s ( ∞ ) = 0, again independent of τ = 0. On the otherhand, we can show that values s ( λ ) of prescribed and random eigenvalues depend on τ and go to0 as τ goes to 0. For this, we need the following lemma. Lemma 5.2.
Let A − λB be a singular pencil and let α be different from all eigenvaluesof 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 ∈ N r ( α ) canbe written as x = q + αq + · · · + α p q p +1 , where the vectors q , . . . , q p +1 form the chain Aq = 0 , Aq = Bq , . . . , Aq p +1 = Bq p , Bq p +1 = 0for certain p ≥
0. Similarly, y ∈ N l ( α ) can be written as y = w + αw + · · · + α r w r +1 , where thevectors w , . . . , w r +1 form the chain A ∗ w = 0 , A ∗ w = B ∗ w , . . . , A ∗ w r +1 = B ∗ w r , B ∗ w r +1 = 0for certain r ≥
0. To show y ∗ Bx = 0 it is enough to show that w ∗ i Bq j = 0 for all i = 1 , . . . , r + 1and j = 1 , . . . , p + 1. For i = 1 or j = p + 1 this follows from w ∗ A = 0 and Bq p +1 = 0, so we canassume that i ≥ j ≤ p . It follows that w ∗ i Bq j = w ∗ i Aq j +1 = w ∗ i − Bq j +1 . As we continuein this manner, we eventually reach either w ∗ A = 0 or Bq p +1 = 0. It follows that y ∗ Bx = 0 andfrom y ∗ Ax = αy ∗ Bx we get y ∗ Ax = 0 as well. Lemma 5.3.
Let λ i be a prescribed or random eigenvalue of (4.1) under the assumptions ofTheorem , where we assume in addition that all prescribed eigenvalues are algebraically simple.Then there exists a positive constant c i such that | s ( λ i ) | = c i | τ | .Proof . First, let λ i be a prescribed eigenvalue. Then by the proof of Theorem 4.6 thecorresponding left and right eigenvectors satisfy y ∗ i ( A − λ i B ) = 0 and ( A − λ i B ) x i = 0 which byLemma 5.2 implies y ∗ i Bx i = 0. But then we have | s ( λ i ) | = c i | τ | with c i = | y ∗ i U D B V ∗ x i | and c i must be nonzero, because λ i is a simple eigenvalue of e A − λ e B for τ = 0.Next, let λ i be a random eigenvalue, such that V ∗ x i = 0 and U ∗ y i = 0. We know (see theproof of Theorem 4.6) that y i is a linear function of τ as y i ( τ ) = w i + τ z i , where w ∗ i ( A − λ i B ) = 0.Since ( A − λ i B ) x i = 0, it follows from Lemma 5.2 that w ∗ i Bx i = 0 and y ∗ i e Bx i = y ∗ i Bx i = τ z ∗ i Bx i .The case V ∗ x i = 0 and U ∗ y i = 0 can be shown analogously.So, if we take a τ of small absolute value and if all finite true eigenvalues are simple and noneof them is too ill-conditioned, then we can separate the finite true eigenvalues from the remainingones using the values s ( λ ).Let ε be the machine precision and let the matrices A and B be scaled in such way that k A k = k B k = 1. If all finite true eigenvalues are simple and not too ill-conditioned, then we OLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION c > τ , andpossibly different for each eigenvalue and each entry in the table. Table 5.1
Characteristics of the eigenvalues of the perturbed pencil as in (4.1) . Eigenvalue λ | s ( λ ) | k V ∗ x k k U ∗ y k Finite true eigenvalue of A − λB c < c ε/ | τ | < c ε/ | τ | Infinite true eigenvalue of A − λB < c ε < c ε/ | τ | < c ε/ | τ | Prescribed eigenvalue of D A − λD B c | τ | c c Random eigenvalue from an L p block c | τ | < c ε/ | τ | c Random eigenvalue from an L Tp block c | τ | c < c ε/ | τ | We now explain the values in Table 5.1. We will start with column | s ( λ ) | and a finite trueeigenvalue, where we assume that all finite true eigenvalues are simple and well-conditioned. Itfollows that λ is a simple eigenvalue of e A − λ e B , therefore y ∗ e Bx = 0 and, since this value isindependent of τ and ε , we have | s ( λ ) | = c . For an infinite eigenvalue we should have y ∗ e Bx = 0in exact computation, instead, in finite precision, we get | y ∗ e Bx | < c ε . Finally, in the genericcase, if λ i is a prescribed or random eigenvalue then Lemma 5.3 yields that | s ( λ i ) | = c i | τ | for apositive constant c i .Finally, the values in the columns k V ∗ x k and k U ∗ x k that are marked by < c ε/ | τ | should bezero in exact arithmetic. In finite precision however, due to the supposed backward stability ofthe applied eigenproblem solver, the computed eigenvalues and eigenvectors of e A − λ e B are exacteigenpairs of a perturbed pencil e A + E − λ ( e B + F ), where k E k ≤ e c k e A k ε and k F k ≤ e c k e B k ε .If we assume that all finite eigenvalues of e A − λ e B are simple, then we have the following resulton the first-order eigenvector perturbations. The proof is omitted since it is a straightforwardgeneralization of the result for the pencil A − λI from [42, Sec. 2.10]. Lemma 5.4.
Let all finite eigenvalues λ i of e A − λ e B be simple and let x i and y i be correspondingleft and right normalized eigenvectors. If the pencil is perturbed into e A + θE − λ ( e B + θF ) , thenthe eigenvector x i perturbs into e x i = x i + θ n X k =1 ,k = i y ∗ k ( E − λ i F ) x i ( λ i − λ k ) s ( λ k ) x k + O ( θ ) . Let λ i be a finite true eigenvalue of A − λB . Then λ i is also an eigenvalue of e A − λ e B and V ∗ x i = 0, where x i is an exact normalized right eigenvector. In finite precision, x i becomesperturbed in the directions of other eigenvectors and by Lemma 5.4 a contribution in the directionof another eigenvector depends on the condition number of the corresponding eigenvalue. Theonly contributions that affect the value of k V ∗ e x i k are those related to prescribed eigenvaluesor random eigenvalues from left singular blocks, as right eigenvectors of other eigenvalues areorthogonal to V . As condition numbers of these eigenvalues are equal to 1 / ( c | τ | ) and k V ∗ x j k = c for the corresponding right eigenvectors, it follows from Lemma 5.4 and the backward stabilityof the computed eigenpairs that k V ∗ e x i k < c ε/ | τ | .Next, we discuss appropriate choices for the value τ . If | τ | is close to ε , then the prescribedand random eigenvalues are very ill-conditioned, and perturbations of eigenvectors may move thevalues of k V ∗ x i k and k U ∗ y i k far away from zero when they should be close to zero. Therefore, if | τ | is too small, we may not be able to use the values of k V ∗ x i k and k U ∗ y i k to extract the trueeigenvalues. Still, if all finite true eigenvalues are simple, then we may use the values | s ( λ i ) | toextract the finite true eigenvalues.8 HOCHSTENBACH, MEHL, AND PLESTENJAK
On the other hand, if | τ | is large, then all eigenvalues, except the infinite ones, are expected tobe well-conditioned which means that the eigenvectors will not change much and the computedleft and right eigenvectors will be orthogonal to V or U in finite precision, when they should be.Therefore, for large | τ | , we can first use max( k V ∗ x i k , k U ∗ y i k ) to extract the true eigenvalues andthen use | s ( λ i ) | to distinguish the finite true eigenvalues from the infinite one. In practice, wesee this as a better option, because it does not depend on finite true eigenvalues being simple.However, we should not choose | τ | too large as this may decrease the precision of the computedfinite true eigenvalues. Since the computed eigenvalues are, due to assumed backward stability,exact eigenvalues of a slightly perturbed pencil e A − λ e B , it is safe to use | τ | up to k e A k ≈ k A k and k e B k ≈ k B k . Also, since from our analysis it follows that only the absolute value of τ seems tomatter, we suggest to choose τ real and positive.Based on the above discussion, we summarize our method in Algorithm 1. Note that wescale the matrices in such way that k A k = k B k = 1, mainly for convenience, to determine anappropriate default value for τ . Algorithm 1: Computing finite eigenvalues of a singular pencil ( A, B ) by a rank-completing perturbation.
Input: A and B , perturbation constant τ (default 10 − ), thresholds δ (default ε / ) and δ (default 10 ε ). Output:
Eigenvalues of the finite regular part. Scale A = (1 /α ) A and B = (1 /β ) B , where α = k A k and β = k B k . Compute nrank(
A, B ): k = rank( A − ζB ) for random ζ . Select random n × k matrices U and V with orthonormal columns. Select diagonal k × k matrices D A and D B such that theeigenvalues of ( D A , D B ) are (likely) different from those of ( A, B )(default: choose diagonal elements of D A and D B uniformly randomfrom the interval [1 , Compute the eigenvalues λ i , i = 1 , . . . , n , and right and lefteigenvectors x i and y i of ( e A, e B ) = ( A + τ U D A V ∗ , B + τ U D B V ∗ ) . Compute s i = y ∗ i e Bx i for i = 1 , . . . , n . Compute ζ i = max( k V ∗ x i k , k U ∗ y i k ) for i = 1 , . . . , n . Return all eigenvalues ( α/β ) λ i , i = 1 , . . . , n , where ζ i < δ and | s i | > δ .As we will show by experiments in the next section, the above approach seems to workvery well in double precision for small or moderate singular pencils. Of course, if some of theeigenvalues are very ill-conditioned (for instance when some of the eigenvalues are multiple), thenthe method may fail in extracting some of the finite true eigenvalues. However, its advantageover staircase-based methods may be the following observation: if we make a wrong rank decisionin a staircase algorithm, then the method usually fails completely and returns no eigenvalues atall; see Example 6.4 in the next section. In contrast, the method proposed here is able to detect,if not all, then at least the well-conditioned finite eigenvalues of the pencil under consideration.
6. Numerical examples.
In this section we demonstrate the method with several numer-ical examples computed in Matlab 2015b. All numerical examples and implementations of the
OLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION
Example 6.1.
We take a 7 × A − λB , where A = − − − − − − −
11 0 0 0 0 0 01 2 1 1 1 1 11 2 3 3 3 3 31 2 3 2 2 2 21 2 3 4 3 3 31 2 3 4 5 5 4 , B = − − − − − − − − − − − − −
12 5 5 5 5 5 52 5 5 4 4 4 42 5 5 6 5 5 52 5 5 6 7 7 72 5 5 6 7 6 6 . The matrices are built in such way that the KCF of the pencil contains blocks of all four possibletypes, nrank(
A, B ) = 6 and the pencil is singular. Its KCF has blocks J (1 / J (1 / N , L ,and L T . If we apply Algorithm 1, we get the values in the following table. Note that values λ k in the first column are values from Line 8 of Algorithm 1, which are scaled back to match theeigenvalues of the original matrix pencil A − λB whose matrices are scaled in Line 1 since theydo not satisfy k A k = k B k = 1. k λ k | s k | k V ∗ x k k k U ∗ y k k . · − . · − . · − . · − . · − . · − ∞ . · − . · − . · − − . . i . · − . · − . · − − . − . i . · − . · − . · − . · − . · − . · − . · − . · − . · − The values in the table follow the pattern from the previous section and it is easy to detectthat λ and λ are finite true eigenvalues, λ is a true infinite eigenvalue, λ , λ , and λ arerandom eigenvalues, and λ is the prescribed eigenvalue. Example 6.2.
We take example C3 from [12] that comes from control theory and belongs toa set of examples C1, C2, and C3, where each has successively more ill-conditioned eigenvalues.The pencil has the form A − λB = − − −
750 0 0 0 2 − λ . Its KCF contains blocks L , J (1), and J (2). As the pencil is rectangular, we add a zero line tomake it square. This adds an L T block to the KCF. Algorithm 1 returns the following table for A − λB , from which the finite true eigenvalues 1 and 2 can be extracted. k λ k | s k | k V ∗ x k k k U ∗ y k k . · − . · − . · − . · − . · − . · − − . . i . · − . · − . · − − . − . i . · − . · − . · − . · − . · − . · − As in [12] we add some noise and perturb initial A − λB into b A − λ b B by adding 10 − rand (4 , A and B . True eigenvalues of b A − λ b B can still be extracted by Algorithm 1 if we adjust theparameter δ . The values we get are in the following table: k λ k | s k | k V ∗ x k k k U ∗ y k k . · − . · − . · − . · − . · − . · − . · − . · − . · − − . . · − . · − . · − . · − . · − . · − HOCHSTENBACH, MEHL, AND PLESTENJAK
Example 6.3.
This is an example from [14, Sec. 5], where the staircase algorithm fails tofind a regular subspace of proper size under a small random perturbation. We take A − λB = − λ δ δ , where δ = 1 . · − . The KCF structure of the pencil is J (0) and L which means that 0 is adouble eigenvalue. It is reported in [14] that if we add a random perturbation of size 10 − to thepencil, then Guptri reports the regular part J (0) and we have been able to confirm this using aMatlab implementation of Guptri in [34]. If we enlarge the perturbation to 10 − , Guptri returnsno regular part at all, while Algorithm 1 returns two finite true eigenvalues λ and λ from thefollowing table. k λ k | s k | k V ∗ x k k k U ∗ y k k − . · − . · − . · − . · − . · − . · − . . · − − . · . · − . · − . · − . · . · − . · − . · − Example 6.4.
We take the singular pencil ∆ − λ ∆ of size 300 ×
300 from [33, Ex. 18].This example is related to two random matrices A and B of size 10 ×
10 in a way that the trueeigenvalues of ∆ − λ ∆ are 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 ∆ − λ ∆ contains 100 N and 10 left and 10 right singular blocks. Theconjecture from [33] is that the singular blocks are 5 L T , 5 L T , 5 L , and 5 L blocks.This example is also available as demo_double_eig_mp in toolbox MultiParEig [35]. Thestaircase algorithm in MultiParEig fails to extract the finite regular part of size 90 in doubleprecision, but manages to extract all 90 finite true eigenvalues using quadruple precision and theMultiprecision Computing Toolbox [36]. If we apply Algorithm 1 to ∆ − λ ∆ in double precision,we get the following values: k λ k | s k | k V ∗ x k k k U ∗ y k k . . i . · − . · − . · − ... ... ... ... ...89 4 . − . i . · − . · − . · − − . . · − . · − . · − ∞ . · − . · − . · − ... ... ... ... ...190 ∞ . · − . · − . · − − . . · − . · − . · − ... ... ... ... ...300 7.125982 2 . · − . · − . · − From the columns k V ∗ x k k and k U ∗ y k k we get max k =1 ,..., (max( k V ∗ x k k , k U ∗ y k k )) = 1 . · − and min k =191 ,..., max( k V ∗ x k k , k U ∗ y k k ) = 4 . · − , which shows a clear gap which separatestrue eigenvalues from the prescribed and random ones. Next, in the set of true eigenvalues thereis also a clear gap between s and s which separates finite true eigenvalues from infinite ones,since min k =1 ,..., | s k | = 3 . · − and max k =91 ,..., | s k | = 1 . · − .
7. The singular two-parameter eigenvalue problem.
We now expand on Section 3.3.In a two-parameter eigenvalue problem (2EP) [1] we have the equations( A + λB + µC ) x = 0 , (7.1) ( A + λB + µC ) x = 0 , OLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION A , B , and C are of size n × n , and A , B , and C are of size n × n . Sought arescalars λ, µ and nonzero vectors x and x such that (7.1) is satisfied. We say that ( λ, µ ) is aneigenvalue of the 2EP and the tensor product x ⊗ x is the corresponding eigenvector. Definethe operator determinants∆ = B ⊗ C − C ⊗ B , ∆ = C ⊗ A − A ⊗ C , (7.2) ∆ = A ⊗ B − B ⊗ A . Then problem (7.1) is related to a coupled pair of GEPs∆ z = λ ∆ z, (7.3) ∆ z = µ ∆ z for a decomposable tensor z = x ⊗ x . If ∆ is nonsingular, then Atkinson [1] shows thatthe solutions of (7.1) and (7.3) agree and the matrices ∆ − ∆ and ∆ − ∆ commute. In thenonsingular case the 2EP (7.1) has n n eigenvalues and it can be solved with a variant of theQZ algorithm on (7.3); see [18].It turns out that for many problems occurring in practice both pencils (∆ , ∆ ) and (∆ , ∆ )are singular and we have a singular 2EP [31]. Applications include delay-differential equations[20], quadratic two-parameter eigenvalue problems [32, 19], model updating [5], and roots ofsystems of bivariate polynomials [37, 2].The eigenvalues of a singular 2EP (7.1) are the finite regular eigenvalues of (7.3); see Sec-tion 3.3. There exists a staircase type algorithm that works on both singular pencils (7.3) simulta-neously and extracts finite regular eigenvalues; see [32] and an implementation in [35]. However,as illustrated in Examples 6.3 and 6.4, a staircase algorithm may fail. In this section we proposean alternative method that may be applied to a singular 2EP, which in some cases finds finiteregular eigenvalues when the staircase algorithm fails, while in some other cases the situation isexactly the opposite.We can apply Algorithm 1 to ∆ z = λ ∆ z , one of the two singular pencils in (7.3), to computethe λ i components of eigenvalues ( λ i , µ i ). This is, however, only half of the required informationand for each λ i we have to find the corresponding µ i . Subsequently, we insert λ = λ i into (7.1) andsearch for common eigenvalues µ of a pair of pencils ( A − λ i B ) − µC and ( A − λ i B ) − µC thatmay be singular as well. We detect the common eigenvalues by comparing the sets of computedeigenvalues for the first and the second pencil, for which we use Algorithm 1 again. The overallmethod is given in Algorithm 2. Algorithm 2: Computing finite regular eigenvalues of a singular 2EPInput:
Matrices A , B , C , A , B , C from (7.1) which provide ∆ and ∆ from (7.2); threshold δ (default δ = ε / ), and parameters for Algorithm 1. Output:
Finite regular eigenvalues of (7.1). Compute finite eigenvalues λ , . . . , λ r of ∆ − λ ∆ using Algo. 1. for j = 1 , . . . , r Compute eigenvalues µ (1)1 , . . . , µ (1) m of ( A − λ j B ) − µC using Algo. 1. Compute eigenvalues µ (2)1 , . . . , µ (2) m of ( A − λ j B ) − µC using Algo. 1. Reorder eigenpairs: | µ (1)1 − µ (2)1 | ≤ · · · ≤ | µ (1) m − µ (2) m | for m = min( m , m ). for k = 1 , . . . , m if | µ (1) k − µ (2) k | < δ then add ( λ j , ( µ (1) k + µ (2) k )) to list of eigenvalues.2 HOCHSTENBACH, MEHL, AND PLESTENJAK
Some remarks about Algorithm 2 are in order. • If we know that each eigenvalue has a unique λ component, then we can replace Lines 6and 7 by selecting ( λ j , ( µ (1)1 + µ (2)1 )) regardless of the difference | µ (1)1 − µ (2)1 | . • If n = n = n then the complexity of Line 1 is O ( n ) while the complexity of Lines 2 to7 is at most O ( n ) in case r = O ( n ). Example 7.1.
Consider a system of bivariate polynomials (cf. [37, Exs. 5.3, 6.2, 6.4]) p ( λ, µ ) = 1 + 2 λ + 3 λ + 4 λ + 5 λµ + 6 µ + 7 λ + 8 λ µ + 9 λµ + 10 µ = 0 ,p ( λ, µ ) = 10 + 9 λ + 8 µ + 7 λ + 6 λµ + 5 µ + 4 λ + 3 λ µ + 2 λµ + µ = 0 . Using a uniform determinantal representation from [2], we write the above system as a 2EP ofthe form A + λB + µC = λ λ − λ
16 + 9 λ + 10 µ − λ − µ − µ ,A + λB + µC = λ λ − λ
15 + 2 λ + µ − λ − µ − µ , where p i ( λ, µ ) = det( A i + λB i + µC i ) for i = 1 ,
2. The obtained 2EP is singular and has 9 regulareigenvalues ( λ 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 computefirst components λ , . . . , λ as finite eigenvalues of the corresponding singular pencil ∆ − λ ∆ from (7.3), whose KCF contains 4 L , 4 L T , 2 N , 1 N , 2 N , and 9 J blocks. For each λ j wecompute the candidates for µ j in Lines 4 and 5, where the KCF of singular pencils ( A i − λ j B i ) − µC i contains 1 N and 3 J blocks for i = 1 , j = 1 , . . . , − λ ∆ might be soill-conditioned that the algorithm cannot separate them from the infinite eigenvalues. In sucha case a possible solution would be to apply computation in higher precision, using, e.g., theMultiprecision Computing Toolbox [36].
8. Conclusions.
We have proposed a method to approximate the finite eigenvalues of asingular pencil by means of a rank-completing perturbation . The use of such a perturbationensures that, generically, the finite and infinite eigenvalues remain fixed, while there appearnewly generated eigenvalues. For many problems we can well distinguish the original eigenvaluesfrom the newly created ones by considering the angles of the eigenvectors with respect to theperturbation spaces, and at the condition numbers of the eigenvalues. Thus, this method may beuseful for a wide range of applications.The proposed method could be an alternative to the class 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, however, the key part of staircase techniques area number of rank decisions, which can be difficult and ill-posed, see e.g., [14] and Examples 6.3and 6.4. In some cases, when these methods fail to return even a single eigenvalue, the newlyproposed method may still compute all or at least some of the eigenvalues.A code for the approach developed in this paper is available in [35].
OLVING SINGULAR GEPs BY A RANK-COMPLETING PERTURBATION Acknowledgments:
The authors would like to thank Stefan Johansson for providing a betaversion of Matrix Canonical Structure (MCS) Toolbox [34] which includes a Matlab implementa-tion of Guptri and is an important alternative to the original Guptri [17] that we can no longeruse in Matlab due to the 32-bit limitation. Furthermore, the authors would like to warmly thanktwo anonymous referees for their careful reading and many expert suggestions and comments ona previous version of this paper.
Genealogical acknowledgment:
During this research project the first two authors found outthat they are twelfth cousins. Christian and Michiel thank their common ancestors CasparH¨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. Plestenjak , Uniformdeterminantal 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 LinearDifferential-Algebraic Equations . In: A. Ilchmann, T. Reis (eds), Surveys in Differential-Elgebraic Equa-tions 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 signalprocessing, 15 (2001), pp. 649–665.[6]
E.J. Davison and S.H. Wang , Properties and calculation of transmission zeros of linear multivariablesystems , Automatica, 10 (1974), pp. 643–658.[7]
F. De Ter´an and F. M. Dopico , Low rank perturbation of Kronecker structures without full rank , SIAMJ. 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 singularmatrix 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: apractical guide , pp. 28–36, SIAM, Philadelphia, 2000.[11]
J. Demmel and B. K˚agstr¨om , Computing stable eigendecompositions of matrix pencils , Linear AlgebraAppl., 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 boundsand 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. Matrix Anal. Appl.,21 (2000), pp. 1004–1025.[15]
A. Emami-Naeini and P. Van Dooren , Computation of zeros of linear multivariable systems,
Automatica18 (1982), pp. 415–430.[16]
F.R. Gantmacher , Theory of Matrices. Volumes 1 and 2 , Chelsea, New York, 1959.[17]
Guptri , software for singular pencils, .[18]
M. E. Hochstenbach, T. Koˇsir, and B. Plestenjak , A Jacobi–Davidson type method for the nonsingulartwo-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-parametereigenvalue problem , Linear Algebra Appl., 436 (2012), pp. 2725–2743.[20]
E. Jarlebring and M. E. Hochstenbach , Polynomial two-parameter eigenvalue problems and matrix pencilmethods 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 eigenvalueof 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, andH. 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 HOCHSTENBACH, MEHL, AND PLESTENJAKPublishing House, Z¨urich, Switzerland, 2006.[24]
A.J. Laub and B.C. Moore , Calculation of transmission zeros using QZ techniques , Automatica, 14 (1978),pp. 557–566.[25] D. S. Mackey, N. Mackey, C. Mehl, and V. Mehrmann , Skew-symmetric matrix polynomials and theirSmith 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. Ma-trix Anal. Appl., 38 (2017), pp. 72–95.[29] ,
Linear algebra properties of dissipative Hamiltonian descriptor systems , SIAM J. Matrix Anal. Appl.,39 (2018), pp. 1489–1519.[30]
V. Mehrmann , private communication, 2018.[31]
A. Muhiˇc and B. Plestenjak , On the singular two-parameter eigenvalue problem , Electron. J. LinearAlgebra, 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, .[35]
MultiParEig . Toolbox for multiparameter eigenvalue problems, .[36]
Multiprecision Computing Toolbox . Advanpix, Tokyo. .[37]
B. Plestenjak and M. E. Hochstenbach , Roots of bivariate polynomial systems via determinantal repre-sentations , 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.), MatrixPencils. Lecture Notes in Mathematics, vol. 973. Springer, Berlin, Heidelberg, 1983, pp. 58–73.[42]
J.H. Wilkinson , The Algebraic Eigenvalue Problem , Oxford University Press, Oxford, 1965.[43] ,
Kronecker’s canonical form and the QZ algorithmalgorithm