
==== Front
BIT Numer Math
BIT Numer Math
Bit. Numerical Mathematics
0006-3835
1572-9125
Springer Netherlands Dordrecht

39022760
1033
10.1007/s10543-024-01033-w
Article
Analysis of eigenvalue condition numbers for a class of randomized numerical methods for singular matrix pencils
Kressner Daniel 1
http://orcid.org/0000-0002-0463-7223
Plestenjak Bor bor.plestenjak@fmf.uni-lj.si

2
1 grid.5333.6 0000000121839049 Institute of Mathematics, EPFL, 1015 Lausanne, Switzerland
2 grid.8954.0 0000 0001 0721 6013 Faculty of Mathematics and Physics, University of Ljubljana, and Institute of Mathematics, Physics and Mechanics, Jadranska 19, 1000 Ljubljana, Slovenia
Communicated by Gunnar Martinsson.

15 7 2024
15 7 2024
2024
64 3 3223 5 2023
4 7 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by/4.0/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
The numerical solution of the generalized eigenvalue problem for a singular matrix pencil is challenging due to the discontinuity of its eigenvalues. Classically, such problems are addressed by first extracting the regular part through the staircase form and then applying a standard solver, such as the QZ algorithm, to that regular part. Recently, several novel approaches have been proposed to transform the singular pencil into a regular pencil by relatively simple randomized modifications. In this work, we analyze three such methods by Hochstenbach, Mehl, and Plestenjak that modify, project, or augment the pencil using random matrices. All three methods rely on the normal rank and do not alter the finite eigenvalues of the original pencil. We show that the eigenvalue condition numbers of the transformed pencils are unlikely to be much larger than the δ-weak eigenvalue condition numbers, introduced by Lotz and Noferini, of the original pencil. This not only indicates favorable numerical stability but also reconfirms that these condition numbers are a reliable criterion for detecting simple finite eigenvalues. We also provide evidence that, from a numerical stability perspective, the use of complex instead of real random matrices is preferable even for real singular matrix pencils and real eigenvalues. As a side result, we provide sharp left tail bounds for a product of two independent random variables distributed with the generalized beta distribution of the first kind or Kumaraswamy distribution.

Keywords

Singular pencil
Singular generalized eigenvalue problem
Eigenvalue condition number
Randomized numerical method
Random matrices
Mathematics Subject Classification

65F15
15A18
15A22
15A21
47A55
68W20
15B52
Swiss National Science Foundation00021 L_192049 Kressner Daniel http://dx.doi.org/10.13039/501100004329 Javna Agencija za Raziskovalno Dejavnost RS N1-0154 P1-0294 Plestenjak Bor issue-copyright-statement© BIT Foundation 2024
==== Body
pmcIntroduction

The purpose of this work is to study three recent numerical methods, introduced in [12, 13], for computing finite eigenvalues of a square singular matrix pencil A-λB, that is, A,B∈Cn×n and det(A-λB)≡0. We say that λ0∈C is an eigenvalue of A-λB if rank(A-λ0B)<nrank(A,B), wherenrank(A,B):=maxζ∈Crank(A-ζB)<n

is the normal rank of the pencil. Similarly, if rank(B)<nrank(A,B) then we say that A-λB has infinite eigenvalue(s).

A major difficulty when working with a singular pencil numerically is the discontinuity of its eigenvalues, that is, the existence of (arbitrarily small) perturbations of A, B that completely destroy the eigenvalue accuracy. To circumvent this phenomenon, it is common to first extract the regular part from the staircase form of the pencil [28] before applying the QZ algorithm [15, 22] to compute the eigenvalues. Notably, the popular software GUPTRI [2, 3] is based on this approach. Numerically, the computation of the staircase form requires several rank decisions and these decisions tend to become increasingly difficult as the algorithm proceeds, which can ultimately lead to a failure of correctly identifying and extracting the regular part [7, 23].

Despite the discontinuity of the eigenvalues mentioned above, Wilkinson [30] observed that the QZ algorithm directly applied to the original singular pencil usually returns the eigenvalues of the regular part with reasonable accuracy. De Terán, Dopico, and Moro [1] explained this phenomenon by developing a perturbation theory for singular pencils, implying that the set of perturbation directions causing discontinuous eigenvalue changes has measure zero. Later on, Lotz and Noferini [18] turned this theory into quantitative statements, by measuring the set of perturbations leading to large eigenvalue changes and defining the notion of weak eigenvalue condition number for singular pencils.

It is important to note that Wilkinson’s observation does not immediately lead to a practical algorithm, because the (approximate) eigenvalues returned by the QZ algorithm are mixed with the spurious eigenvalues originating from the (perturbed) singular part and it is a nontrivial task to distinguish these two sets. During the last few years, several methods have been proposed to circumvent this difficulty.

Inspired by the findings in [1], Hochstenbach, Mehl, and Plestenjak [12] proposed to introduce a modification of the form1 A~-λB~:=A-λB+τ(UDAV∗-λUDBV∗)

for matrices DA,DB∈Ck×k with k:=n-nrank(A,B), random matrices U,V∈Cn×k, and a scalar τ≠0. Generically, A~-λB~ is a regular pencil and the regular part of A-λB is exactly preserved, i.e., if λi is a finite eigenvalue of A-λB then λi is an eigenvalue of A~-λB~ with the same partial multiplicities. More specifically, λi is an eigenvalue of A-λB if and only if λi is an eigenvalue of A~-λB~ such that its right/left eigenvectors x, y satisfy V∗x=0 and U∗y=0. The latter property is used to extract the eigenvalues of A-λB from the computed eigenvalues of A~-λB~.

In [13], two different variations of the approach from [12] described above are proposed. Instead of adding a modification, the pencil is projected to the generically regular pencil U⊥∗AV⊥-λU⊥∗BV⊥ for random matrices U⊥,V⊥∈Cn×(n-k), and the eigenvalues of A-λB are extracted from the computed eigenvalues of the smaller pencil. The third method analyzed in this work consists of computing the eigenvalues of A-λB from the augmented generically regular pencilAUTASAV∗0-λBUTBSBV∗0,

where SA,SB,TA,TB∈Ck×k and U,V∈Cn×k are random matrices. For both variants, it holds generically that the regular part of A-λB is fully preserved in exact arithmetic. Note that roundoff error affects this eigenvalue preservation property when forming the modified pencils (1) and U⊥∗AV⊥-λU⊥∗BV⊥ in finite precision.

One goal of this work is to show that the modifications introduced by the three methods above are numerically safe. More specifically, we show that, with high probability, the eigenvalue condition numbers of the modified pencils are not much larger than the weak eigenvalue condition numbers of the original pencil. In particular, these methods can be expected to return good accuracy for well-conditioned eigenvalues of A-λB in the presence of roundoff error. Another implication of our result is that the eigenvalue condition numbers of the modified pencils represent a reliable complementary criterion for identifying reasonably well-conditioned finite eigenvalues in any of the algorithms from [12, 13].

Related work Closer to the analyses in [1, 18], it was recently suggested in [17] to perturb the full pencil: A+τE-λ(B+τF), where E,F∈Cn×n are random Gaussian matrices and τ>0 is small but well above the level of machine precision. Unlike for the three methods mentioned above, the regular part of A-λB is not preserved by this perturbation. On the other hand, the direct connection to [18] allows to facilitate their analysis and use the computed eigenvalue condition numbers of the perturbed pencil as a criterion to identify finite eigenvalues of the original pencil. In [17, P. 2], it was stated that a similar analysis would be more difficult for the method from [12] because of the structure imposed on the random perturbation in (1). In this work, we will address this question and carry over the analysis from [17, 18] to the three methods above. In particular, our analysis confirms that the computed eigenvalue condition numbers can be used as a reliable indicator for such algorithms as well.

Outline The structure of the paper is as follows. In Sect. 2 we review basic concepts for singular pencils as well as δ-weak condition numbers. In Sect. 3 we present the three randomized numerical methods that we analyze in Sect. 4, where we also obtain the new left tail bounds. This is followed by numerical examples in Sect. 5. In the appendix we provide results obtained with symbolic computation that verify the results from latter sections.

Preliminaries

Reducing subspaces and eigenvectors

In order to define eigenvectors of a singular pencil according to [1, 18], we first introduce the Kronecker canonical form (KCF) and the notion of minimal reducing subspaces, see, e.g., [9, 29].

Theorem 1

(Kronecker canonical form) Let A,B∈Cn×n. Then there exist nonsingular matrices P,Q∈Cn×n such that2 P(A-λB)Q=R(λ)00S(λ),R(λ)=J-λIr00Is-λN,

where J and N are in Jordan canonical form with N nilpotent. Furthermore,S(λ)=diag(Lm1(λ),⋯,Lmk(λ),Ln1(λ)T,⋯,Lnk(λ)T),

where Lj(λ)=[0Ij]-λ[Ij0] is of size j×(j+1), and mi,ni≥0 for i=1,⋯,k, where k=n-nrank(A,B).

The pencil R(λ) in (2) is called the regular part of A-λB and contains the eigenvalues of A-λB, where the Jordan part J-λIr contains the Jordan blocks of the finite eigenvalues of A-λB. If λ0 is a finite eigenvalue of A-λB, then its partial multiplicities are the sizes of the Jordan blocks associated with λ0. A finite eigenvalue is called simple if it is a simple root of detR(λ). The pencil S(λ) is called the singular part of A-λB and contains right singular blocks Lm1(λ),⋯,Lmk(λ) and left singular blocks Ln1(λ)T,⋯,Lnk(λ)T, where m1,⋯,mk and n1,⋯,nk are called the right and left minimal indices of the pencil, respectively.

We say that a subspace M is a reducing subspace [29] for the pencil A-λB if dim(AM+BM)=dim(M)-k, where k=n-nrank(A,B) counts the number of right singular blocks. The minimal reducing subspace MRS(A,B) is the intersection of all reducing subspaces and is spanned by the columns of Q corresponding to the blocks Lm1(λ),…,Lmk(λ). Analogously, L is a left reducing subspace for the pencil A-λB if dim(A∗L+B∗L)=dim(L)-k and the minimal left reducing subspace LRS(A,B) is the intersection of all left reducing subspaces.

For an eigenvalue λ0∈C of A-λB, a nonzero vector x∈Cn is called a right eigenvector if (A-λ0B)x=0 and x∉MRS(A,B). A nonzero vector y∈Cn such that y∗(A-λ0B)=0 and y∉LRS(A,B) is called a left eigenvector. This agrees with the definition of eigenvectors from [5, 18]. Compared to a regular pencil, the eigenvectors of singular pencils have a much larger degree of non-uniqueness, due to components from the minimal reducing subspaces.

Eigenvalue perturbation theory

Suppose that we perturb an n×n matrix pencil A-λB into3 A~-λB~:=A+ϵE-λ(B+ϵF),

where ϵ>0. We define ‖(E,F)‖:=(‖E‖F2+‖F‖F2)1/2, where ‖·‖F denotes the Frobenius norm of a matrix. When ‖(E,F)‖=1 we can identify the pencil E-λF with a point on the unit sphere in C2n2 and think of (E, F) as a direction of the perturbation (3).

Before addressing the singular case, let us first recall the classical eigenvalue perturbation theory [27] for a regular pencil A-λB. Consider a simple finite eigenvalue λ0∈C of A-λB with normalized right/left eigenvectors x, y. For ‖(E,F)‖=1 and sufficiently small ϵ>0 there exists an eigenvalue λ0(ϵ) of the perturbed pencil (3) satisfying the perturbation expansion4 |λ0(ϵ)-λ0|=|y∗(E-λ0F)x||y∗Bx|ϵ+O(ϵ2)≤(1+|λ0|2)1/2|y∗Bx|ϵ+O(ϵ2)

as ϵ→0. Note that the inequality becomes an equality for the direction E=(1+|λ0|2)-1/2yx∗, F=-λ¯0E. In turn, the absolute condition number of λ0, defined as5 κ(λ0)=limϵ→0sup‖(E,F)‖≤11ϵ|λ0(ϵ)-λ0|,

satisfies6 κ(λ0)=1/γ(λ0),whereγ(λ0)=|y∗Bx|(1+|λ0|2)-1/2;

see, e.g., [8, Lemma 3.1 and Eq. (3.3)].

For a singular pencil, the definition (5) always leads to an infinite condition number because of the discontinuity of eigenvalues. To address this, we first recall the eigenvalue expansion by De Terán, Dopico, and Moro [1, Corollary 2].

Theorem 2

Let λ0 be a finite simple eigenvalue of an n×n pencil A-λB of normal rank n-k, k≥1. Let X=[X1x], Y=[Y1y] be n×(k+1) matrices with orthonormal columns such that: X1 is a basis for ker(A-λ0B)∩MRS(A,B), X is a basis for ker(A-λ0B), Y1 is a basis for ker((A-λ0B)∗)∩LRS(A,B), and Y is a basis for ker((A-λ0B)∗). If E-λF is such that det(Y1∗(E-λ0F)X1)≠0, then, for sufficiently small ϵ>0, there exists an eigenvalue λ0(ϵ) of the perturbed pencil (3) such that7 λ0(ϵ)=λ0-det(Y∗(E-λ0F)X)y∗Bx·det(Y1∗(E-λ0F)X1)ϵ+O(ϵ2).

The expansion (7) allows one to define the directional sensitivity for perturbations E-λF, ‖(E,F)‖=1, satisfying the condition of Theorem 2:8 σE,F(λ0)=det(Y∗(E-λ0F)X)y∗Bx·det(Y1∗(E-λ0F)X1);

see [18, Definition 2.4 and Corollary 3.3]. We can now generalize the definition of γ(λ0) in (6) to singular pencils by using the right/left eigenvectors x, y from Theorem 2. Because of the appearance of the factor y∗Bx in the denominator of (8), we expect that the quantity 1/γ(λ0) continues to play a crucial role in determining the sensitivity of an eigenvalue. However, it is important to note that x, y are particular choices of eigenvectors made in Theorem 2: the right eigenvector x is orthogonal to the right minimal reducing subspace MRS and the left eigenvector y is orthogonal to LMR. Under these constraints, the eigenvectors x and y of the simple eigenvalue λ0 become uniquely determined up to multiplication by unit complex numbers. It follows from Y1∗BX1=0, Y1∗Bx=0 and y∗BX1=0 that this particular choice maximizes |y∗Bx|, i.e., for all vectors z and w such that ‖z‖2=‖w‖2=1, (A-λ0B)z=0, and w∗(A-λ0B)=0, it holds that9 |w∗Bz|(1+|λ0|2)-1/2≤γ(λ0).

Remark 3

If A, B are real matrices and λ0 is a real simple eigenvalue, then the matrices X,X1,Y,Y1 and vectors x, y in Theorem 2 can be chosen to be real as well.

Clearly, the perturbations for which the quantity det(Y1∗(E-λ0F)X1) in the denominator of (7) vanishes form a set of measure zero in the unit sphere in C2n2. In other words, this event has zero probability if we draw (E, F) uniformly at random from the the unit sphere in C2n2, which we will denote by (E,F)∼U(2n2). The δ-weak condition number introduced by Lotz and Noferini [18, Definition 2.5] offers a more refined picture by measuring the tightest upper bound t such that the directional sensitivity (8) stays below t with probability at least 1-δ.

Definition 1

Let λ0∈C be a finite simple eigenvalue of a singular pencil A-λB. The δ-weak condition number of λ0 is defined asκw(λ0;δ)=inf{t∈R:P(σE,F(λ0)<t)≥1-δ},(E,F)∼U(2n2).

If λ0 is a finite simple eigenvalue of a regular pencil A-λB, it follows from (3) and (6) that σE,F(λ0)=|y∗(E-λ0F)x|/|y∗Bx|≤κ(λ0)=1/γ(λ0) for all ‖(E,F)‖=1. Therefore, if we apply Definition 1 to a regular pencil, it follows from (5) that κw(λ0;δ) converges to κ(λ0)=1/γ(λ0) monotonically from below as δ↓0. The following result from [18, Theorem 5.1] and [17, Theorem 3.1] suggests to use 1/γ(λ0) as a proxy for eigenvalue sensitivity in the singular case as well.

Theorem 4

Let λ0∈C be a finite simple eigenvalue of an n×n singular pencil A-λB of normal rank n-k. Then for δ≤k/(2n2) it holds that1δ2n2·1γ(λ0)≤κw(λ0;δ)≤kδ2n2·1γ(λ0),

where γ(λ0) is defined as in (6) with the right/left eigenvectors x, y from Theorem 2.

The algorithms for the singular generalized eigenvalue problem presented in the next section use γ(λ0) to identify (finite) eigenvalues numerically.

We finish this section by remarking that the absolute condition number (5) is consistent with the definitions in [17, 18]. Following, e.g., [10] one could also consider a notion of relative condition number that imposes ‖E‖F≤‖A‖F and ‖F‖F≤‖B‖F instead of ‖(E,F)‖≤1 on the perturbation direction in (3). In effect, both the standard and weak (absolute) condition numbers get multiplied by the factor (‖A‖F2+‖B‖F2)1/2. Under reasonable choices of the parameters involved, this additional factor does not differ significantly for the modified pencils. In particular, our results presented for absolute condition numbers easily extend to relative condition numbers.

Randomized numerical methods based on the normal rank

In this section, we describe in some detail the three numerical methods from [12, 13] for computing the finite eigenvalues of an n×n singular pencil A-λB with nrank(A,B)=n-k for k≥1. All three methods require knowledge about the exact normal rank in order to leave the regular part intact. If this quantity is not known a priori, it can be determined from rank(A-ξiB) for a small number of randomly chosen ξi∈C.1

Rank-completing modification

We first consider the rank-completing method from [12], where a random pencil of normal rank k is added to yield a (generically regular) matrix pencil10 A~-λB~:=A-λB+τ(UDAV∗-λUDBV∗),

where DA,DB∈Ck×k are diagonal matrices such that DA-λDB is regular, U,V∈Cn×k are matrices of rank k, and τ∈C is nonzero. Note that k is the smallest normal rank for such a modification to turn a singular into a regular pencil. The following result [12, Summary 4.7], see also [13, Remark 3.5], characterizes the dependence of eigenvalues and eigenvectors of the modified pencil (10) on τ, DA, DB, U, and V∗.

Summary 5

Let A-λB be an n×n singular pencil of normal rank n-k with left minimal indices n1,⋯,nk and right minimal indices m1,⋯,mk. Let N=n1+⋯+nk and M=m1+⋯+mk. Then the regular part of A-λB has size r:=n-N-M-k and generically (with respect to the entries of DA,DB,U,V∗), the modified pencil A~-λB~ defined in (10) is regular and its eigenvalues are classified in the following four disjoint groups: True eigenvalues: There are r eigenvalues counted with their multiplicities that coincide with the eigenvalues of the original pencil A-λB with the same partial multiplicities. The right/left eigenvectors x, y of A~-λB~ belonging to these eigenvalues satisfy the orthogonality relations V∗x=0 and U∗y=0.

Prescribed eigenvalues: There are k eigenvalues such that V∗x≠0 for all right eigenvectors x and U∗y≠0 for all left eigenvectors y. These are the k eigenvalues of DA-λDB.

Random right eigenvalues: There are M eigenvalues, which are all simple and such that V∗x=0 for all right eigenvectors x and U∗y≠0 for all left eigenvectors y.

Random left eigenvalues: There are N eigenvalues, which are all simple and such that V∗x≠0 for all right eigenvectors x and U∗y=0 for all left eigenvectors y.

Summary 5 has the following practical consequences. If we compute all eigenvalues λi of (10), together with the (normalized) right and left eigenvectors xi and yi for i=1,…,n, then max(‖V∗xi‖2,‖U∗yi‖2)=0 if and only if λi is an eigenvalue of A-λB. In numerical computations, we can use max(‖V∗xi‖2,‖U∗yi‖2)<δ1, where δ1 is a prescribed threshold, as a criterion to extract the true eigenvalues in the first phase. Note that for a simple finite eigenvalue xi and yi are unique (up to multiplication by unit complex numbers) because, generically, the modified pencil is regular. They correspond to eigenvectors of the original singular pencil satisfying the orthogonality constraints V∗xi=0 and U∗yi=0.

In the second phase, we use the (reciprocal) eigenvalue sensitivities for extracting simple finite eigenvalues, that is, we compute11 γi=|yi∗Bxi|(1+|λi|2)-1/2

and identify λi as a finite eigenvalue if γi>δ2 for a prescribed threshold δ2. Note that 1/γi is the absolute condition number of λi as an eigenvalue of the (generically) regular pencil (10); see (6).

For different matrices U and V in (10) we obtain different eigenvectors xi and yi and thus different values of γi for the same eigenvalue, while the changing of τ, DA and DB does not affect the eigenvectors [13, Lemma 3.4]. In Sect. 4 we will analyze these values for random U and V and compare them to the unique value γ(λi) that appears in the δ-weak condition number of λi as an eigenvalue of the singular pencil A-λB; see Theorem 4.

The considerations above lead to Algorithm 1 from [12]. In theory, the results returned by the algorithm are independent of τ≠0. In practice, |τ| should be neither too small nor too large in order to limit the impact of roundoff error; in [12] it is suggested to scale A and B so that ‖A‖1=‖B‖1=1, where ‖A‖1=max1≤j≤n∑i=1n|aij|, and take τ=10-2. The quantity ε stands for the machine precision.

By the theory in [12, 13], the eigenvectors of the modified pencil (10) that correspond to true eigenvalues do not change if we replace U and V by U~=UR and V~=VS, where R and S are arbitrary nonsingular k×k matrices. Thus, choosing U, V to have orthonormal columns does not violate the genericity assumption in Summary 5.

Note that γi in line 4 was initially computed in [12] as γi=|yi∗B~xi|. This was changed to γi=|yi∗B~xi|(1+|λi|2)-1/2 in [13] to be consistent with [17, 18]. Since for true eigenvalues yi∗B~xi=yi∗Bxi, we use B instead of B~ to simplify the analysis.

We remark that the values γi from (11) were also used in [17] for computing finite eigenvalues of a singular pencil via unstructured random perturbations. The use of full-rank perturbations comes with two disadvantages: the orthogonality relations for the eigenvectors exploited above are not satisfied and, in contrast to Algorithm 1, the eigenvalues of the perturbed pencil in [17] differ from the exact eigenvalues of the original pencil. The latter leads one to choose τ>0 very small, but at the same time it needs to stay well above the level of machine precision.

Normal rank projections

In [13], a variant of Algorithm 1 was proposed that uses random projections to a pencil of smaller size, equal to the normal rank n-k. In addition, the method does not require to choose the matrices DA, DB, and the parameter τ.

For U⊥,V⊥∈Cn×(n-k), we consider the (n-k)×(n-k) pencil U⊥∗(A-λB)V⊥. To connect it to the modified pencil (10) used in Algorithm 1, let us assume that the columns of U⊥ and V⊥ span the orthogonal complements of ranges of U and V, respectively, so that U⊥∗U=0, V⊥∗V=0. The pencilA^-λB^:=[UU⊥]∗(A-λB+τ(UDAV∗-λUDBV∗)[VV⊥]=U∗(A-λB)V+τ(DA-λDB)U∗(A-λB)V⊥U⊥∗(A-λB)VU⊥∗(A-λB)V⊥

is then equivalent to (10) and we observe the following.

Proposition 1

([13, Proposition 4.1 and Theorem 4.2]) Let A-λB be a complex n×n singular pencil of normal rank n-k. Then, under the assumptions of Summary 5, the (n-k)×(n-k) pencil A22-λB22:=U⊥∗(A-λB)V⊥ is generically regular and the eigenvalues of A22-λB22 are precisely: the random eigenvalues of (10) (groups 3 and 4 in Summary 5);

the true eigenvalues of A-λB.

Based on the above results, an algorithm is devised in [13]. Algorithm 2 is a simplified form that matches Algorithm 1 as much as possible.

The following corollary shows that the reciprocal eigenvalue condition number γi computed in line 4 of Algorithm 2 matches the corresponding quantity of Algorithm 1.

Corollary 1

Let λi∈C be a simple eigenvalue of a singular pencil A-λB. Under the assumptions of Proposition 1, if (λi,xi,yi) is an eigentriple of (10) such that ‖xi‖2=1, ‖yi‖2=1, V∗xi=0, and U∗yi=0; and

(λi,wi,zi) is an eigentriple of A22-λB22 such that ‖wi‖2=1, ‖zi‖2=1, U∗(A-λiB)V⊥wi=0, and zi∗U⊥∗(A-λiB)V=0,

then |yi∗Bxi|=|zi∗U⊥∗BV⊥wi|.

Proof

Since λi is simple, the vectors xi,yi from (a) and wi,zi from (b) are uniquely defined up to multiplication by unit complex numbers. If (a) and (b) both hold then it immediately follows that, up to sign changes, xi=V⊥wi and yi=U⊥zi. □

Augmentation

The third method, also presented in [13], uses the (n+k)×(n+k) augmented (or bordered) matrix pencil12 Aa-λBa:=AUTASAV∗0-λBUTBSBV∗0,

where SA,SB,TA, and TB are k×k diagonal matrices and U, V are n×k matrices.

Proposition 2

([13, Proposition 5.1]) Let A-λB be an n×n singular pencil of normal rank n-k such that all its eigenvalues are semisimple. Assume that the diagonal k×k pencils SA-λSB and TA-λTB are regular and that their 2k eigenvalues are pairwise distinct. Furthermore, let U,V∈Cn×k have orthonormal columns such that the augmented pencil (12) is regular. Then the pencil (12) has the following eigenvalues: 2k prescribed eigenvalues, which are precisely the eigenvalues of SA-λSB and TA-λTB;

the random eigenvalues of (10) (groups 3 and 4 in Summary 5) with the same U and V and with DA=TASA, DB=TBSB;

the true eigenvalues of A-λB.

The algorithm based on the above proposition is given in Algorithm 3.

Again, the reciprocal eigenvalue condition number γi computed in line 5 of Algorithm 3 matches the corresponding quantity of Algorithm 1.

Corollary 2

Under the assumptions of Proposition 2, if (λi,xi,yi) is an eigentriple of (10) such that ‖xi‖2=1, ‖yi‖2=1, V∗xi=0, and U∗yi=0; and

(λi,wi,zi), where wi=[wi1T0]T and zi=[zi1T0]T, such that wi1∈Cn, zi1∈Cn, ‖wi‖2=1, ‖zi‖2=1, is an eigentriple of the augmented pencil (12),

then |yi∗Bxi|=|zi∗Bawi|=|zi1∗Bwi1|.

Proof

If a) and b) are both true then it immediately follows that, up to sign changes, xi=wi1 and yi=zi1.□

Let us remark that Algorithms 1–3, with additional heuristic criteria, can in practice, due to the “positive” effect of roundoff error, successfully compute multiple finite eigenvalues as well; we refer to [13, Sec. 6] for details.

Remark 6

In exact arithmetic, a simple finite eigenvalue λ0 of a singular pencil A-λB is an exact eigenvalue of a modified regular pencil A~-λB~ that we obtain by rank-completion, projection, or augmentation. In finite precision, A~-λB~ is affected by roundoff error. First, the construction of A~-λB~ itself already introduces some error:13 (A~+▵A~)-λ(B~+▵B~).

Standard error analysis results [11] imply for all three methods that ‖▵A~‖2≲u‖A‖2 and ‖▵B~‖2≲u‖B‖2, where u denotes unit roundoff, provided that the parameters are reasonably chosen. For example, in the augmentation (12), the norms of TA,SA should be balanced with the norm of A and the norms of TB,SB should be balanced with the norm of B. The application of the QZ algorithm for computing the eigenvalues of the perturbed pencil (13) introduces further error but, thanks to its backward stability, the additional error is of the same nature. In summary, the computed eigenvalues are the exact eigenvalues of a slightly perturbed modified pencil. In turn, we expect the computed approximation for an eigenvalue to be accurate if the (exact) eigenvalue of A~-λB~ is not ill-conditioned. Our analysis will show that if λ0 is well-conditioned as an eigenvalue of A-λB, then it is – with high probability – also a well-conditioned eigenvalue of A~-λB~. In turn, one can expect good accuracy for well-conditioned eigenvalues of A-λB.

Let us stress that it is nontrival to relate (13) back to an error in the original singular pencil A-λB. At least, we are not aware of existing backward error analysis results for singular matrix pencils that would allow us to do this.

Probabilistic analysis

Our goal is to analyze the behavior of the quantities γi in Algorithms 1–3 and show that they are unlikely to be much below γ(λi). It follows from Corollaries 1 and 2 that it is sufficient to consider Algorithm 1 and the quantity γi defined in (11).

In the following, we assume that λi is a simple eigenvalue of A-λB and let X=[X1x] and Y=[Y1y] denote the orthonormal bases for ker(A-λiB) and ker((A-λiB)∗) introduced in Theorem 2. We recall from Theorem 4 that the reciprocal of γ(λi)=|y∗Bx|(1+|λi|2)-1/2 critically determines the sensitivity of λi as an eigenvalue of A-λB. The sensitivity of λi as an eigenvalue of A~-λB~ from (10) is given by 1/γi with γi=|yi∗Bxi|(1+|λi|2)-1/2; see (11). The eigenvectors xi,yi are normalized (‖xi‖2=‖yi‖2=1) and depend on the choices of U and V in Algorithm 1. Generically (with respect to DA,DB,U,V∗), Summary 5.1 yields the relations14 xi=[X1x]aαandyi=[Y1y]bβ,V∗xi=0,U∗yi=0.

If V∗X1 is invertible, the choice of V entirely determines |α| because V∗xi=0 implies a=-(V∗X1)-1V∗xα and, hence, ‖xi‖2=1 implies|α|=1/1+‖(V∗X1)-1V∗x‖22.

Analogously, the choice of U determines |β| if U∗Y1 is invertible.

Since Y∗BX1=0 and Y1∗BX=0, we get yi∗Bxi=αβ¯y∗Bx and thusγi=|α||β|γ(λi).

The relation 0≤|α|,|β|≤1 immediately gives γi≤γ(λi), in line with (9). A small value of |α||β| means an increased eigenvalue sensitivity for A~-λB~, potentially causing Algorithm 1 to yield unnecessarily inaccurate results. In the following, we will show that this is unlikely when random matrices U, V are used in Algorithm 1. This also implies that the (reciprocal) condition numbers computed in Algorithm 1 can be used with high probability to correctly identify finite simple eigenvalues.

Preliminary results

Let N1(μ,σ2) denote the normal distribution with mean μ and variance σ2. In particular, x∼N1(0,1) is a standard (real) normal random variable. We write z∼N2(0,1) if z=x+iy is a standard complex normal variable, that is, x,y∼N1(0,12) are independent. In the following, we will analyze real matrices (F=R) and complex matrices (F=C) simultaneously. For this purpose, we set ϕ=1 for F=R and ϕ=2 for F=C.

The matrices U and V from Algorithm 1 belong to the Stiefel manifoldVkn(F)={Q∈Fn×k:Q∗Q=I}.

We will choose them randomly (and independently) from the uniform distribution on Vkn(F). A common way to compute such a matrix is to perform the QR decomposition of an n×k Gaussian random matrix M, that is, the entries of M are i.i.d. real or complex standard normal variables, see e.g., [21, 26]. That this indeed yields the uniform distribution follows from the following variant of the well-known Bartlett decomposition theorem ([24, Theorem 3.2.14], [6, Proposition 7.2]); see also [18, Proposition 4.5].

Theorem 7

For F∈{R,C}, let M∈Fn×k, n≥k, be a Gaussian random matrix. Consider the QR decomposition M=QR, where Q∈Vkn(F) and R∈Fk×k is upper triangular with non-negative diagonal entries. Then the entries of Q and the entries of the upper triangular part of R are all independent random variables;

Q is distributed uniformly over Vkn(F);

rij∼Nϕ(0,1) for 1≤i<j≤k;

ϕrjj2∼χ2(ϕ(n-j+1)) for j=1,…,k;

where χ2(ℓ) denotes the chi-squared distribution with ℓ degrees of freedom.

Part (b) of Theorem 7 implies that each column of Q is distributed uniformly over the unit sphere in Fn. Note that x=z/‖z‖2, for a Gaussian random vector z∈Fn, has the same distribution. The following result provides the distribution of the entries of x; this result can be found for F=R in [4].

Lemma 1

Consider a random vector x distributed uniformly over the unit sphere in Fn for n≥2. Then the entries of x are i.i.d. with|xi|2∼Betaϕ2,ϕ(n-1)2,i=1,…,n,

where Beta denotes the beta distribution.

Proof

By Theorem 7, the entries of x are independent. Without loss of generality, let i=1. Using that x=z/‖z‖2 for a Gaussian random vector z and setting w=z2…zn, it follows that |x1|2=|z1|2|z1|2+‖w‖22, where z1, w are independent and ϕ|z1|2∼χ2(ϕ), ϕ‖w‖22∼χ2(ϕ(n-1)). This implies the claimed result; see, e.g., [6, p. 320]. □

Our analysis will connect α and β from (14) to the nullspaces of k×(k+1) standard Gaussian matrices, which are characterized by the following result.

Lemma 2

For a Gaussian random matrix Ω∈Fk×(k+1) with k≥2, let x be a vector in the nullspace of Ω such that ‖x‖2=1. Then, with probability one, |xi| is uniquely determined and satisfies15 |xi|2∼Betaϕ2,ϕk2.

Proof

We assume that Ω has rank k, which holds with probability 1. For a Gaussian random vector ω∈Fk+1 independent of Ω∗, consider the QR decomposition [Ω∗,ω]=QR. Letting x denote the last column of Q, it follows that x is orthogonal to the columns of Ω∗ or, in other words, x is in the nullspace of Ω. From Theorem 7, it follows that x is distributed uniformly over the unit sphere in Fk+1. The distribution (15) then follows from Lemma 1. Finally, note that |xi|2 is uniquely determined because the nullspace of Ω has dimension 1. □

Statistics of |α|,|β|

The results above readily yield the distribution of |α|2 and |β|2.

Proposition 3

For F∈{R,C}, let U, V be n×k independent random matrices from the uniform distribution on the Stiefel manifold Vkn(F). Consider a finite simple eigenvalue λi∈F of a singular pencil A-λB with A,B∈Fn×n. Let xi and yi be the right and left normalized eigenvectors of the perturbed regular pencil (10) for the eigenvalue λi and let α,β be defined as in (14). Then |α| and |β| are independent random variables and|α|2,|β|2∼Betaϕ2,ϕk2,

where ϕ=1 for F=R and ϕ=2 for F=C.

Proof

We will only prove the distribution for α; the derivation for β is entirely analogous. By the unitary invariance of the uniform distribution over the Stiefel manifold we may assume without loss of generality that X1=[e1…ek] and x=ek+1, where ei is the i-th vector of the standard canonical basis. By Theorem 7, the matrix V is obtained from the QR factorization Ω=VR of an n×k Gaussian random matrix Ω, with R being invertible almost surely. We partitionV∗=V1v2⋯=R-∗Ω1ω2⋯=R-∗Ω∗,

such that V1,Ω1 are k×k matrices and v2,ω2 are vectors. Then0=V∗xi=V∗X1a+V∗xα=V1v2aα=R-∗Ω1ω2aα.

Since submatrices of Gaussian random matrices are again Gaussian random matrices, this means that aα is in the nullspace of a k×(k+1) Gaussian random matrix and has norm 1. Thus, the result on the distribution of α follows from Lemma 2.

The independence of |α| and |β| follows from the independence of U and V combined with the fact that |α| does not depend on U and |β| does not depend on V. □

Remark 8

It is important to emphasize that the case F=R in Proposition 3 not only requires A,B,DA,DB to be real but also the eigenvalue λi to be real.

Remark 9

An analysis similar to the one above was performed in [18, Proposition 6.5] and [17, Section 4.1] for unstructured perturbations. This analysis also starts from the relation (14) and then analyzes the distribution of |α|·|β|. One significant difference in our case is that α and β are independent due to the structure of the perturbation in (10), while this does not hold for the setting considered in [17, 18].

Statistics of |α|·|β|

As explained above, we aim at showing that the random variable |α||β| is unlikely to become tiny. We start by computing the expected value of |α||β|. Since |α| and |β| are independent random variables, we have E[|α||β|]=E[|α|]E[|β|]. The factors can be computed using the following result from [18, Lemma A.1].

Lemma 3

Let X∼Beta(a,b), where a,b>0. ThenE[X1/2]=B(a+1/2,b)B(a,b),whereB(a,b)=Γ(a)Γ(b)Γ(a+b).

To simplify the presentation, we will from now on denote the scalars α and β, in the setting of Proposition 3, as αF and βF with F∈{R,C}. Combining Proposition 3 and Lemma 3 gives the following result.

Lemma 4

Under the assumptions of Proposition 3, the following holds: E[|αC|]=E[|βC|]=πΓ(k+1)2Γ(k+3/2) and E|αC||βC|=πΓ(k+1)24Γ(k+3/2)2.

E[|αR|]=E[|βR|]=Γ((k+1)/2)πΓ((k+2)/2) and E|αR||βR|=Γ((k+1)/2)2πΓ((k+2)/2)2.

Table 1 contains the computed expected values for different k, using the results of Lemma 4.

Using the well-known bounds16 x≤Γ(x+1)Γ(x+1/2)≤x+1/2

for x>0, the expected values of |α||β| from Lemma 4 can be bounded as17 π4(k+1)≤E[|αC||βC|]≤π4(k+1/2),2π(k+1)≤E[|αR||βR|]≤2πk.

Table 1 Expected values of |α||β| for k=n-nrank(A,B)=1,2,4,…,64

k	E|αC||βC|	E|αR||βR|	
1	0.44444	0.40528	
2	0.28444	0.25000	
4	0.16512	0.14063	
8	0.08972	0.07477	
16	0.04688	0.03857	
32	0.02398	0.01959	
64	0.01213	0.00987	

One therefore expects that γi underestimates the true values γ(λi) by roughly a factor 1/k.

For real matrices A and B, one would prefer to use real matrices in the perturbation (10) as well, because eigenvalue computations are performed more efficiently in real arithmetic. As we can see from Table 1 as well as from the bounds (17), the expected value of |α||β| for real perturbations is only slightly smaller than the one for complex perturbations. However, as we will see in the following, the left tail of |α||β| is less favorable in the real case and it appears to be safer to use complex modifications of the original pencil even for real data.

Bounds on left tail of |α|·|β|

We will start with a simple tail bound that extends a result from [17, Proposition 4.3] to the complex case.

Corollary 3

Under the assumptions of Proposition 3, we have18 P(|αC||βC|<t)≤2kt

and19 P(|αR||βR|<t)≤8kt/π

for every 0≤t≤1.

Proof

Using min{|α|4,|β|4}≤|α|2|β|2, we obtainP|α|2|β|2<t2≤Pmin|α|4,|β|4<t2≤P|α|2<t+P|β|2<t=2B(ϕ/2,ϕk/2)∫0txϕ/2-1(1-x)ϕk/2-1dx,

where ϕ=1 for F=R and ϕ=2 for F=C. For ϕ=2, we obtain from |1-x|≤1 thatP|αC|2|βC|2<t2≤2tB(1,k)=2kt.

Similarly, for ϕ=1 we getP|αR|2|βR|2<t2≤4tB(1/2,k/2)≤8ktπ,

where we used the bound (16) to derive the last inequality. □

Although we know from Proposition 3 that |α| and |β| are independent random variables, this is not used in Corollary 3. The following proposition, where we exploit the fact that |α| and |β| are independent, significantly improves the results of Corollary 3.

Proposition 4

Under the assumptions of Proposition 3, it holds that20 P|αC||βC|<t≤k2t2(1-2lnt)

and21 P|αR||βR|<t≤2kπt(-lnt)+O(t),t→0.

Proof

To derive the bound (20) we first observe that, since |αC| and |βC| are independent by Proposition 3, it holds thatP(|αC||βC|≤t)=P|αC|2|βC|2≤t2=∬Dg(x;k)g(y;k)dxdy

with D=(x,y)∈[0,1]×[0,1]:xy≤t2 and g(x;k)=1B(1,k)(1-x)k-1. Using the bound g(x;k)≤1B(1,k)=k gives (20):P(|αC||βC|≤t)≤k2∬Ddxdy=k2t2(1-2lnt).

For real perturbations we have |αR|2,|βR|2∼Beta(1/2,k/2) with the distribution function satisfyingh(x;k)=1B(1/2,k/2)x-1/2(1-x)k/2-1≤1B(1/2,k/2)x-1/2.

This givesP(|αR||βR|≤t)=∬Dh(x;k)h(y;k)dxdy≤1B(1/2,k/2)2∬Dx-1/2y-1/2dxdy=1B(1/2,k/2)24t2π-1+t-lnt≤2kπt2π-1+t-lnt,

where we applied (16). □

Proposition 4 indicates that, even for a real singular pencil, it would be better to use complex perturbations as they give a much smaller probability of obtaining tiny |α||β|. This is confirmed by the following lower bound for P(|αR||βR|<t).

Lemma 5

Under the assumptions of Proposition 3, it holds that22 P|αR||βR|<t≥8(k-1)πt+O(t2),t→0.

Proof

It is easy to see thatP(|αR||βR|<t)≥P(|αR|<t)+P(|βR|<t)-P(|αR|<t)P(|βR|<t)=2P(|αR|<t)-P(|αR|<t)2.

From the Taylor expansion of P(|αR|<t) around t=0 we getP(|αR||βR|<t)≥4B(1/2,k/2)t+O(t2)

and the bound follows from applying (16). □

From Lemma 5 and Proposition 4 we see that for sufficiently small t>0 the lower bound (22) for P|αR||βR|<t is much larger than the upper bound (20) for P|αC||βC|<t.

Figure 1 compares the obtained bounds to P|αC||βC|<t and P|αR||βR|<t, computed for k=4 and k=8 using the probability density functions from Appendix A. The solid black line corresponds to P|αC||βC|<t. The blue dotted and dashed lines are the refined bound (20) from Proposition 4 and the simple upper bound (18) from Corollary 3, respectively. The solid red line corresponds to P|αR||βR|<t. The corresponding bounds are magenta curves, which show the lower bound (22) from Lemma 5, the refined upper bound (21) from Proposition 4, and the simple upper bound (19) from Corollary 3, respectively, As expected, the bounds from Proposition 4 are much sharper than the simple bounds from Corollary 3. Also, it is clearly seen from Fig. 1 that the probability of obtaining a tiny value for |α||β| is much larger when using real perturbations.Fig. 1 Comparison of bounds and actual values of P(|α||β|<t) for t=10-5,10-4,10-3,10-2, and for k=4 (left) and k=8 (right)

Summary Together with the discussion in the beginning of this section, Proposition 4 allows us to compare the (reciprocal) eigenvalue sensitivity γ(λi) of the original pencil with the corresponding quantity γi for any of the three modified pencils used in Algorithms 1–3. For complex perturbations, we obtain thatγi≤γ(λi)≤γi/t

holds with probability at least 1-k2t2(1-2lnt) for any t>0.

Numerical examples

All numerical examples were obtained with Matlab 2021b [19]. We used the implementations of Algorithms 1–3 available as routine singgep in MultiParEig [25].

Example 1

We consider the 8×8 singular matrix pencil23 A=-1-1-1-1-1-1-10100000001211111012333330123222201234333-11234554100002212,B=-2-2-2-2-2-2-202-1-1-1-1-1-1025555550255444402556555-1255677712556766100000-1-10,

which is constructed so that the KCF contains blocks of all four possible types. It holds that nrank(A,B)=6, and the KCF has blocks J1(1/2), J1(1/3), N1, L0, L1, L0T, and L2T.

Algorithm 1 was applied 105 times using random real and complex modifications; we compared the computed values of γ1 to the exact value of γ(λ1) for the eigenvalue λ1=1/3. The histograms of γ1/γ(λ1) for real and complex modifications together with the corresponding probability density function (pdf) from Appendix A are presented in Fig. 2. The histograms appear to be consistent with the pdfs. The computed average values of |α||β| are 0.28437 for complex and 0.24934 for real modifications, which are both close to the theoretically predicted values for k=2 in Table 1. We note that we get almost identical results for the other eigenvalue λ2=1/2.

Fig. 2 Example 1: Histogram of γ1/γ(λ1) and pdf using complex (left) and real (right) modifications

Example 2

For the second example we consider the pencil Δ1-λΔ0 from [12, Ex. 7.1] with matricesΔ0=B1⊗C2-C1⊗B2,Δ1=C1⊗A2-A1⊗C2

of size 25×25 related to the two-parameter eigenvalue problem (for details see, e.g., [12, Sec. 7]) of the formA1+λB1+μC1=004+7λ1005+8λ2-λ16+9λ+10μ310-λ1-μ00001-μ00,A2+λB2+μC2=007+4λ1006+3λ9-λ15+2λ+μ8100-λ1-μ00001-μ00,

wheredet(A1+λB1+μC1)=1+2λ+3μ+4λ2+5λμ+6μ2+7λ3+8λ2μ+9λμ2+10μ3,det(A2+λB2+μC2)=10+9λ+8μ+7λ2+6λμ+5μ2+4λ3+3λ2μ+2λμ2+μ3.

The normal rank of Δ1-λΔ0 is 21 and the KCF contains 4 L0, 4 L0T, 2 N4, 1 N2, 2 N1, and 9 J1 blocks. Its finite eigenvalues are λ-components of the 9 solutions of the system of two bivariate polynomials det(A1+λB1+μC1)=0 and det(A2+λB2+μC2)=0. The pencil Δ1-λΔ0 has one real eigenvalue λ1=-2.41828 and eight complex eigenvalues λ2,…,λ9.

Algorithm 2 was applied 105 times using random real and complex projections. We compared the computed values of γ1 to the exact value of γ(λ1) for the real eigenvalue λ1. In both cases the method successfully computed all finite eigenvalues. The histograms of γ1/γ(λ1) together with the corresponding pdf from Appendix A are presented in Fig. 3. Similar to the previous example, both histograms appear to be consistent with the pdfs. The computed average values of |α||β| are 0.16413 for complex and 0.14124 for real projections; both are again very close to the values in Table 1 for k=4.

Fig. 3 Example 2: Histogram of γ1/γ(λ1) and pdf using complex (left) and real (right) projections

Let us note that we get essentially identical results if we exchange the methods and use Algorithms 2 or 3 in Example 1 and Algorithms 1 or 3 in Example 2, as expected by Corollary 1. However, if we apply Algorithm 1 or 3 with real modifications or Algorithm 2 with real projections to the real pencil Δ1-λΔ0, and consider any of the 8 complex eigenvalues λ2,…,λ9, we get results that cannot be explained by Proposition 3 and behave like the results in the next example.

Example 3

In this example we study the effect of real projections on a complex eigenvalue of a real pencil, which is a situation that is not covered by Proposition 3. By considering two equivalent singular pencils we will show that the distribution function for |αR||βR|, when we apply real projections to complex eigenvalues of real pencils, depends on more than just the difference k between the size of the pencil and its normal rank, which is the key value in Proposition 3.

We take block diagonal matrices24 A0=e1e5TA20,B0=e1e4TB20,

where e1,e4,e5 are standard basis vectors in R5, andA20=10000020000011000-11000001,B20=1000001000001000001000000.

The 10×10 pencil A0-λB0 has nrank(A0,B0)=6 and the KCF has blocks J1(1+i), J1(1-i), J1(2), N1, 3 L0, L1, 3 L0T, and L1T. We multiply A0 and B0 into Ai=QiA0Zi and Bi=QiB0Zi by real matrices Qi and Zi whose entries are independent random variables uniformly distributed on (0, 1) for i=1,2 to get two equivalent singular pencils of the same size, normal rank and eigenvalues.

Algorithm 2 was applied 105 times using random real projections to each of the pencils A1-λB1 and A2-λB2 and the computed value of γ1 was compared to the exact value of γ(λ1) for the complex eigenvalue λ1=1+i. The histograms of γ1/γ(λ1) together with the theoretical distribution functions from Sect. A are presented in Figure 4. We see that, although the pencils A1-λB1 and A2-λB2 are equivalent, the histograms are different. The histograms also look different than in the case when we use real projections for a real eigenvalue for k=4 (Fig. 2 right) or complex projections (Fig. 2 left). While the shape of the left histogram in Fig. 4 resembles the shape expected for complex perturbations, is the shape of the right histogram more in line with the distribution function for real perturbations. The computed average values of |αR||βR| are also different, we get 0.12195 for A1-λB1 and 0.13623 for A2-λB2, both values are completely different from the values in Table 1 for k=4.

Since we get histograms of different shape for two real equivalent pencils, this shows that the distribution function for |αR||βR|, when we apply real perturbations to complex eigenvalues of real pencils, depends on more than just the structure of the KCF.Fig. 4 Example 3: Histogram of γ1/γ(λ1) and theoretical distribution functions (real and complex) for k=4 using real perturbations for a complex eigenvalue of a real pencil for A1-λB1 (left) and A2-λB2 (right)

We remark that the histograms (not shown) for the real eigenvalue λ3=2 for both pencils look identical to the right picture in Fig. 3 from Example 2, where k=4 as well, which agrees with Remark 3 that even if some eigenvalues of the real pencil are complex, this does not affect the behavior of real perturbations to real eigenvalues. The computed average values of |αR||βR| for the real eigenvalue λ3 are 0.14077 for A1-λB1 and 0.14110 for A2-λB2, they both agree with the value in Table 1 for k=4.

The last numerical example reflects that the case of real perturbations for a complex eigenvalue of a real singular pencil is not covered by Proposition 3. Indeed, this case was not taken properly into account in [18] and it remains an open problem to derive an expression or a tight simple bound for the δ-weak condition number of a complex eigenvalue under real perturbations.

Conclusions

We have analyzed three random based numerical methods for computing finite eigenvalues of a singular matrix pencil. All algorithms are based on random matrices that transform the original singular pencil into a regular one in such way that the eigenvalues remain intact. Our analysis confirms the numerical validity of these methods with high probability.

We also obtained sharp left tail bounds on the distribution of a product of two independent random variables distributed with the generalized beta distribution of the first kind or Kumaraswamy distribution.

A closer look at the distribution of |α||β|

In this section, we derive explicit expressions for the probability density function (pdf) of |α||β|. For this purpose, we first write down the pdfs of |α|, |β|, which we can express with the generalized beta distribution of the first kind [20] in the real case and the Kumaraswamy distribution [14] in the complex case.

Lemma 6

Under the assumptions of Proposition 3, we have |αC|,|βC|∼GB1(2,1,1,k)=Kumaraswamy(2,k) with the pdf h(x;k)=2kx(1-x2)k-1,

|αR|,|βR|∼GB1(2,1,12,k2) with the pdf g(x;k)=2B(12,k2)(1-x2)k/2-1,

where GB1(a,b,p,q) is the generalized beta distribution of the first kind and Kumaraswamy(p,q) denotes the Kumaraswamy distribution with parameters p and q.

Proof

The pdf of X∼Beta(p,q) is given byf(x;p,q)=1B(p,q)xp-1(1-x)q-1.

It follows that the pdf for Y=X1/2 is2yf(y2;p,q)=2B(p,q)y2p-1(1-y2)q-1.

Now, we apply Proposition 3 and insert the corresponding values of p, q. The proof is concluded by identifying the obtained densities with the Kumaraswamy distribution in the complex case and the generalized beta distribution of the first kind in the real case; see, e.g., [14, 20]. □

The refs [14, 20] used in the proof above also provide expressions for the moments of the two distributions in Lemma 6, which can be used to verify the results of Lemma 4.

In order to derive the pdf of |α|·|β| from the result of Lemma 6, we have performed symbolic computations in Wolfram Mathematica [31]. Let fC(x;k) and fR(x;k) denote the pdfs for |αC||βC| and |αR||βR|, respectively. Using that |αC|,|βC| are independent and |αC|,|βC|∼Kumaraswamy(2,k), we obtain thatfC(x;k)=2k2x(x2-1)2k-1B(k,k+1)(k(x2-1)2k+22F1(k+1,k+1,2k+2,1-x2)-22F1(k,k,2k+1,1-x2)),

where 2F1(a,b,c,z) is the hypergeometric function. One finds that fC takes the form25 fC(x;k)=x(pk-1(x2)+qk-1(x2)lnx),

where pk-1 and qk-1 are polynomials of degree k-1. Some explicit expressions for small k arefC(x;1)=-4xlnx,fC(x;2)=16x(-1+x2-(1+x2)lnx),fC(x;3)=18x(-3+3x4-2(1+4x2+x4)lnx),fC(x;4)=323x(-11-27x2+27x4+11x6-6(1+9x2+9x4+x6)lnx).

We were unable to obtain a closed form for the distribution of |αR||βR|, with |αR|,|βR| independent and |αR|,|βR|∼GB1(2,1,12,k2). For k=2m we conjecture thatfR(x;2m)=pm-1(x2)+qm-1(x2)lnx,

where pm-1 and qm-1 are polynomials of degree m-1. In addition to the fact that this looks similar to (25), this conjecture is also supported by the following expressions for small m:fR(x;2)=-lnx,fR(x;4)=94(-1+x2-(1+x2)lnx),fR(x;6)=225128(-3+3x4-2(1+4x2+x4)lnx).

The graphs of fC(x;k) and fR(x;k) for k=2,4,6,8 are presented in Fig. 5. Note that limx→0fC(x;k)=0 and limx→0fR(x;k)=∞.Fig. 5 Pdfs fC(x;k) (left) and fR(x;k) (right) for k=2,4,6,8

Acknowledgements

We thank Ivana Šain Glibić for discussions related to this work. The authors also thank the reviewers for their careful reading and useful comments.

Funding

Funding was provided by Swiss National Science Foundation (Grant No. 00021 L_192049) and Javna agencija za znanstvenoraziskovalno in inovacijsko dejavnost Republike Slovenije (Grant Nos. N1-0154 and P1-0294).

Declarations

Conflict of interest

Not applicable.

1 Although this heuristics works very well in practice, an analysis is beyond the scope of this work. Already the seemingly simpler special case of numerically deciding whether A-λB is close to a singular pencil is quite intricate [16].

The work of the first author was supported by the SNSF research project Probabilistic methods for joint and singular eigenvalue problems, Grant Number: 200021 L_192049. The work of the second author was supported by the Slovenian Research and Innovation Agency (Grants N1-0154, P1-0294).

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. De Terán F Dopico FM Moro J First order spectral perturbation theory of square singular matrix pencils Linear Algebra Appl. 2008 429 2–3 548 576
De Terán, F., Dopico, F.M., Moro, J.: First order spectral perturbation theory of square singular matrix pencils. Linear Algebra Appl. 429(2–3), 548–576 (2008)
2. Demmel J Kågström B The generalized Schur decomposition of an arbitrary pencil A-λB: robust software with error bounds and applications. I. Theory and algorithms ACM Trans. Math. Softw. 1993 19 2 160 174
Demmel, J., Kågström, B.: The generalized Schur decomposition of an arbitrary pencil : robust software with error bounds and applications. I. Theory and algorithms. ACM Trans. Math. Softw. 19(2), 160–174 (1993)
3. Demmel J Kågström B The generalized Schur decomposition of an arbitrary pencil A-λB: robust software with error bounds and applications. II. Software and applications ACM Trans. Math. Softw. 1993 19 2 175 201
Demmel, J., Kågström, B.: The generalized Schur decomposition of an arbitrary pencil : robust software with error bounds and applications. II. Software and applications. ACM Trans. Math. Softw. 19(2), 175–201 (1993)
4. Dixon JD Estimating extremal eigenvalues and condition numbers of matrices SIAM J. Numer. Anal. 1983 20 4 812 814
Dixon, J.D.: Estimating extremal eigenvalues and condition numbers of matrices. SIAM J. Numer. Anal. 20(4), 812–814 (1983)
5. Dopico FM Noferini V Root polynomials and their role in the theory of matrix polynomials Linear Algebra Appl. 2020 584 37 78
Dopico, F.M., Noferini, V.: Root polynomials and their role in the theory of matrix polynomials. Linear Algebra Appl. 584, 37–78 (2020)
6. Eaton ML Multivariate Statistics 1983 New York Wiley
Eaton, M.L.: Multivariate Statistics. Wiley, New York (1983)
7. Edelman A Ma Y Staircase failures explained by orthogonal versal forms SIAM J. Matrix Anal. Appl. 2000 21 3 1004 1025
Edelman, A., Ma, Y.: Staircase failures explained by orthogonal versal forms. SIAM J. Matrix Anal. Appl. 21(3), 1004–1025 (2000)
8. Frayssé V Toumazou V A note on the normwise perturbation theory for the regular generalized eigenproblem Numer. Linear Algebra Appl. 1998 5 1 1 10
Frayssé, V., Toumazou, V.: A note on the normwise perturbation theory for the regular generalized eigenproblem. Numer. Linear Algebra Appl. 5(1), 1–10 (1998)
9. Gantmacher, F.R.: The Theory of Matrices. Vols. 1, 2. Chelsea Publishing Co., New York (1959). Translated by K. A. Hirsch
10. Higham DJ Higham NJ Structured backward error and condition of generalized eigenvalue problems SIAM J. Matrix Anal. Appl. 1999 20 2 493 512
Higham, D.J., Higham, N.J.: Structured backward error and condition of generalized eigenvalue problems. SIAM J. Matrix Anal. Appl. 20(2), 493–512 (1999)
11. Higham NJ Accuracy and Stability of Numerical Algorithms 2002 2 Philadelphia SIAM
Higham, N.J.: Accuracy and Stability of Numerical Algorithms, 2nd edn. SIAM, Philadelphia (2002)
12. Hochstenbach ME Mehl C Plestenjak B Solving singular generalized eigenvalue problems by a rank-completing perturbation SIAM J. Matrix Anal. Appl. 2019 40 3 1022 1046
Hochstenbach, M.E., Mehl, C., Plestenjak, B.: Solving singular generalized eigenvalue problems by a rank-completing perturbation. SIAM J. Matrix Anal. Appl. 40(3), 1022–1046 (2019)
13. Hochstenbach ME Mehl C Plestenjak B Solving singular generalized eigenvalue problems. Part II: projection and augmentation SIAM J. Matrix Anal. Appl. 2023 44 4 1589 1618
Hochstenbach, M.E., Mehl, C., Plestenjak, B.: Solving singular generalized eigenvalue problems. Part II: projection and augmentation. SIAM J. Matrix Anal. Appl. 44(4), 1589–1618 (2023)
14. Jones MC Kumaraswamy’s distribution: a beta-type distribution with some tractability advantages Stat. Methodol. 2009 6 1 70 81
Jones, M.C.: Kumaraswamy’s distribution: a beta-type distribution with some tractability advantages. Stat. Methodol. 6(1), 70–81 (2009)
15. Kågström B Kressner D Multishift variants of the QZ algorithm with aggressive early deflation SIAM J. Matrix Anal. Appl. 2006 29 1 199 227
Kågström, B., Kressner, D.: Multishift variants of the QZ algorithm with aggressive early deflation. SIAM J. Matrix Anal. Appl. 29(1), 199–227 (2006)
16. Kressner D Voigt M Distance problems for linear dynamical systems Numerical Algebra. Matrix Theory, Differential-Algebraic Equations and Control Theory 2015 Cham Springer 559 583
Kressner, D., Voigt, M.: Distance problems for linear dynamical systems. In: Numerical Algebra. Matrix Theory, Differential-Algebraic Equations and Control Theory, pp. 559–583. Springer, Cham (2015)
17. Kressner D Šain Glibić I Singular quadratic eigenvalue problems: linearization and weak condition numbers BIT 2023 63 1 18, 25
Kressner, D., Šain Glibić, I.: Singular quadratic eigenvalue problems: linearization and weak condition numbers. BIT 63(1), 18, 25 (2023)
18. Lotz M Noferini V Wilkinson’s bus: weak condition numbers, with an application to singular polynomial eigenproblems Found. Comput. Math. 2020 20 6 1439 1473
Lotz, M., Noferini, V.: Wilkinson’s bus: weak condition numbers, with an application to singular polynomial eigenproblems. Found. Comput. Math. 20(6), 1439–1473 (2020)
19. Matlab: The MathWorks Inc., Natick, MA, Version 2021b
20. McDonald JB Some generalized functions for the size distribution of income Econometrica 1984 52 647 664
McDonald, J.B.: Some generalized functions for the size distribution of income. Econometrica 52, 647–664 (1984)
21. Mezzadri F How to generate random matrices from the classical compact groups Notices Am. Math. Soc. 2007 54 5 592 604
Mezzadri, F.: How to generate random matrices from the classical compact groups. Notices Am. Math. Soc. 54(5), 592–604 (2007)
22. Moler CB Stewart GW An algorithm for generalized matrix eigenvalue problems SIAM J. Numer. Anal. 1973 10 241 256
Moler, C.B., Stewart, G.W.: An algorithm for generalized matrix eigenvalue problems. SIAM J. Numer. Anal. 10, 241–256 (1973)
23. Muhič A Plestenjak B A method for computing all values λ such that A+λB has a multiple eigenvalue Linear Algebra Appl. 2014 440 345 359
Muhič, A., Plestenjak, B.: A method for computing all values such that has a multiple eigenvalue. Linear Algebra Appl. 440, 345–359 (2014)
24. Muirhead RJ Aspects of Multivariate Statistical Theory. Wiley Series in Probability and Mathematical Statistics 1982 New York Wiley
Muirhead, R.J.: Aspects of Multivariate Statistical Theory. Wiley Series in Probability and Mathematical Statistics. Wiley, New York (1982)
25. Plestenjak, B.: MultiParEig. Toolbox for multiparameter and singular eigenvalue problems. www.mathworks.com/matlabcentral/fileexchange/47844-multipareig
26. Stewart, G.W.: The efficient generation of random orthogonal matrices with an application to condition estimators. SIAM J. Numer. Anal. 17(3), 403–409 (1980). ((loose microfiche suppl.))
27. Stewart GW Sun J-G Matrix Perturbation Theory 1990 New York Academic Press
Stewart, G.W., Sun, J.-G.: Matrix Perturbation Theory. Academic Press, New York (1990)
28. Van Dooren P The computation of Kronecker’s canonical form of a singular pencil Linear Algebra Appl. 1979 27 103 140
Van Dooren, P.: The computation of Kronecker’s canonical form of a singular pencil. Linear Algebra Appl. 27, 103–140 (1979)
29. Van Dooren P Kågström B Ruhe A Reducing subspaces: Definitions, properties and algorithms Matrix Pencils 1983 New York Springer 58 73
Van Dooren, P.: Reducing subspaces: Definitions, properties and algorithms. In: Kågström, B., Ruhe, A. (eds.) Matrix Pencils, pp. 58–73. Springer, New York (1983)
30. Wilkinson JH Kronecker’s canonical form and the QZ algorithm Linear Algebra Appl. 1979 28 285 303
Wilkinson, J.H.: Kronecker’s canonical form and the algorithm. Linear Algebra Appl. 28, 285–303 (1979)
31. Wolfram Research Inc., Champaign, IL. Mathematica, Version 13.0
